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is  allowed  to  be  time  sequential  it  performs  uniformly  better  than  the  ad-hoc 
detector.  Furthermore,  the  optimum  detector  demonstrates  an  adaptive  learning 
quality,  through  the  process  of  updating  the  a-priori  probability  density  functions 
of  the  unknown  parameters,  from  iteration  to  iteration. 
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Introduction 


Acoustic  reverberation  is  a  phenomenon  which  active  sonar  systems  have 
to  combat.  Reverberation  results  when  the  propagating  pulse’s  energy  is  scattered 
from  inhomogeneities  in  the  ocean  and  its  boundaries.  Reverberation  data  usually 
has  a  complex  and  highly  variable  range-Doppler  map.  which  depends  on  the 
location  of  the  sonar  array  in  the  water  column,  and  on  its  beampatterns. 

Ocean  boundary  reverberation  (surface  or  bottom),  can  leak  through 
sidelobes  or  even  through  the  edges  of  a  main  lobe  in  a  conventional  forward 
looking  beam  output.  The  above  leakage  typically  has  a  sudden  onset  and  can 
contribute  Lj  the  overall  noise  level.  In  conventional  active  sonar  systems, 
discrimination  against  boundary  reverberation  is  typically  achieved  through 
controlling  the  sidelobe  characteristics  of  the  main,  forward  looking  beam. 

A  more  sophisticated  approach  to  the  problem  has  suggested  the  use  of 
an  adaptive  beamforming  structure  which  dynamically  steers  spatial  nulls  in  the 
direction  of  the  interferent  boundary  patch  (surface  or  bottom).  Such  a  structure 
is  called  in  the  literature  a  noise  canceler66.  In  addition  to  the  conventional 
forward  looking  beam,  the  adaptive  noise  canceling  structure  uses  one  or  more 
reference  beams.  The  reference  beams  are  constructed  such  that  they  receive  well 
the  boundary  reverberation,  but  are  essentially  prevented  from  receiving  signals 
coming  in  from  the  main  look  direction.  The  processor  adaptively  filters  the 
reference  beams  to  provide  a  good  estimate  of  the  boundary  reverberation 
contaminant  in  the  main  beam.  It  then  subtracts  the  adaptive  filter  output  from 
the  main  beam  output  to  get  the  desired  -  reverberation  canceled  -  signal. 
Implementations  of  such  adaptive  beamforming  structures  which  take  advantage 
of  the  spatial  separation  between  desired  signals  and  boundary  reverberation 
contaminants  are  presented  in23.  Additional  relevant  sources  are54-64. 
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The  above  solution  to  the  problem  is  representative  of  an  approach  which 
is  based  on  intuition,  but  it  is  not  at  all  clear  that  it  yields  an  optimal  processor 
(under  any  optimization  criterion). 

Another  approach  suggests  treating  the  problem  as  a  whole  right  from  the 
beginning  without  imposing  intuitive  components  on  the  processor  structure,  and 
using  all  a-priori  knowledge  available.  Detection  theory  provides  us  with  a 
mathematical  framework  out  of  which  optimum  processors  can  be  designed.  The 
processor  will  evolve  out  of  the  mathematical  solution  of  the  problem,  and  will 
not  be  restricted  to  using  familiar  structures. 

Although  Bayes  optimal  processors  have  been  derived  for  the  case  of 
volume  reverberation35,  little  work  has  been  done  which  takes  advantage  of  a- 
priori  knowledge  of  the  time-evolving  spatial  characteristics  of  boundary 
reverberation.  Related  Bayes  optimal  work  concerning  interference  sources  of 
certain  and  uncertain  (but  not  time  varying)  location  is  contained  in  1  and  26-24-25 
respectively. 

Here,  a  classical  detection  theoretic  approach  will  be  applied  to  the 
processing  of  a  vector  time  series.  That  vector  may  be  composed  of  the  single 
array  element  outputs,  or  of  some  preformed  beams.  Some  optimality  criterion  is 
chosen,  and  then  the  processor  structure  is  allowed  to  evolve  freely  out  of  the 
mathematical  solution  of  the  problem.  Any  uncertain  parameters  are  treated  as 
random  variables  and  all  knowledge  about  them  is  summarized  in  a-priori 
probability  density  functions. 

This  dissertation  is  organized  as  follows:  Chapter  1  reviews  the  current 
scientific  knowledge  about  the  reverberation  phenomenon.  It  then  details  some  of 
the  physical  and  statistical  properties  of  reverberation  important  for  attempting 
to  quantify  spatial  covariance  matrix  expressions  of  it.  These  matrices  play  a 
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major  role  in  the  optimum  detector  expressions. 

In  the  second  chapter,  an  important  tool  is  developed.  Often  times,  when 
likelihood  ratio  expressions  are  complex,  computing  the  detector  performance 
analytically  is  either  very  hard,  or  impossible.  Then,  being  capable  of  synthesizing 
multi-element  (or  multi-beam)  active  sonar  pings  becomes  a  prerequisite  to  one's 
ability  to  test  detection  algorithms  using  a  Monte-Carlo  simulation  approach.  In 
lieu  of  sonar  synthesis  programs,  which  are  time-consuming  and  impractical  for 
this  application,  this  chapter  uses  a  spectral  estimation  approach,  to  compute  the 
multi-channel,  time  evolving  power  spectrum  of  the  problem,  including  the 
transmit  pulse,  the  element  (beam)  transmit  and  receive  beampatterns,  the  volume 
and  boundary  effects,  and  the  dynamics  of  the  scenario.  Then,  the  multi-channel 
covarianc**  matrices  are  computed,  and  the  normal  equations  problem  is  solved, 
yielding  a  time  variable,  one-step  forward  prediction  filter.  The  inverse  filter 
frequency  response  is  an  estimate  of  the  spectral  character  of  the  problem.  Now, 
all  that's  left  is  to  pass  vectors  of  uncorrelated  white  noise  through  the  filter,  in 
order  to  get  multiple  sonar  pings  of  the  same  sonar  scenario,  for  a  relatively  low 
time  investment.  Lastly,  the  synthesis  capability  is  demonstrated  through  an 
example. 

Chapter  3  reviews  the  two  approaches  to  detection  of  active  sonar  pings 
in  a  reverberating  environment.  As  discussed  above,  the  first  is  a  combination  ad- 
hoc  noise  canceler-matched  filter.  The  chapter  then  lays  the  giound  idles  for 
Bayes  optimum  detection,  and  examines  the  array  detection  problem.  Lastly,  it 
develops  the  Bayes  optimum  detector  (BOD)  for  the  simple  case  of  signal  known 
exactly  in  noise  known  exactly,  which  turns  out  to  be  a  matched  filter. 

All-zero  adaptive  filters  exist  in  many  variations,  where  the  main 
difference  is  in  the  adaptive  algorithm  modifying  the  filter  weights.  In  the  fourth 
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chapter,  the  Joint  Process  Least  Squares  Lattice  (JCLSL)  all-zero  filter  is  reviewed, 
and  then  its  pole-zero  (HR)  counterpart  is  developed  using  an  embedding 
approach,  which  makes  use  of  the  all-zero  filter  recursions  to  create  the  pole-zero 
filter  recrsions.  Then,  through  a  set  of  simulations,  it  is  shown  that  in  some 
scenarios,  it  is  beneficial  for  the  adaptive  filter  to  possess  the  ability  to  adapt 
poles.  Another  simulation  shows  that  in  the  scenarios  used  in  this  dissertation,  the 
performance  advantage  of  the  pole-zero  filter  is  marginal. 

The  fifth  chapter  develops  the  Bayes  optimum  detectors  for  the  three 
cases  tested  in  this  dissertation: 

1.  Signal  known  exactly  in  correlated  Gaussian  interference  of  precisely  known 
direction. 

2.  Signal  known  exactly  in  correlated  Gaussian  interference  of  uncertain 
direction. 

3.  Signal  known  exactiy  in  correlated  Gaussian  interference  of  uncertain 
direction,  plus  a  fixed  interference  of  uncertain  direction. 

The  sixth  chapter  distinguishes  between  two  different  implementations 
for  the  Bayes  optimum  detector,  namely  block  processing  and  time  sequential.  It 
develops  the  time-sequential  detector  design  equations  including  the  updates  of 
the  a-priori  probability  density  function  estimates  which  increase  its  knowledge 
about  the  problem  as  time  progresses. 

Chapter  7  reviews  the  main  methods  used  to  evaluate  and  compare  the 
performance  of  the  various  detection  schemes. 

Chapters  8  compares  the  performance  of  both  detection  approaches,  when 
the  Bayes  optimum  detector  is  restricted  to  the  block  processing  method.  Finally, 
chapter  9  concludes  this  work  by  comparing  the  performance  of  both  detectors 
when  this  time  the  BOD  is  allowed  the  time  sequential  method.  It  also 
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demonstrates  the  adaptive  quality  of  the  time  sequential  BOD  through  the 
updating  process  of  the  a-priori  probability  density  function  from  one  iteration  to 
the  next. 


1.  Reverberation  in  the  Ocean: 

Active  underwater  electronic  systems  propagate  acoustic  signals  through 
the  water  column  and  use  the  information  carried  by  these  signals  to  detect 
targets,  echo  range,  or  communicate  voice  and  data.  Like  every  other  electronic 
system,  these  systems  suffer  from  noise  effects,  which  mask  the  desired  return 
signal,  and  limit  their  performance.  Underwater  systems’  performance  in  general  is 
limited  by  external  noise  sources  and  not  internal  (thermal)  noise.  In  the  absence 
of  hostile  jammers,  this  external  noise  is  usually  divided  into  two  main  categories, 
namely  ambient  ocean  noise,  and  reverberation.  Among  ambient  ocean  noise 
contributors  one  can  find  rain,  waves,  seismic  noise,  various  biological  noise 
sources,  jhip  traffic  etc.,  and  they  ail  exist  independent  of  the  presence  of  the 
underwater  system.  An  inherent  difference  between  reverberation  noise  and 
ocean  noise  is  that  reverberation  exists  only  as  a  result  of  acoustic  energy  being 
transmitted,  whereas  ambient  ocean  noise  is  always  there. 

System’s  performance  may  be  limited  by  either  type  of  interference 
(reverberation  or  sea  noise),  and  in  fact  the  same  system  may  be  reverberation 
limited  for  part  of  its  operating  range,  and  ambient  noise  limited  for  the  rest  of  its 
operating  range.  The  systems  considered  in  this  d.ssertation  are  assumed  to 
operate  in  conditions  which  dictate  that  performance  be  reverberation  limited. 
Therefore,  ambient  ocean  noise  is  neglected  in  all  derivations. 

Reverberation  is  the  result  of  scattering  of  energy  originating  from  the 
propagating  pulse,  by  inhomogeneities  in  the  ocean  and  its  boundaries.  In  some 
respects  this  problem  is  akin  to  the  radar  "clutter"  problem.  Reverberation  is 
usually  divided  into  three  classes,  namely,  surface,  volume  and  bottom 
reverberation.  Irregularities  in  the  ocean  surface  and  the  acoustic  impedance 
contrast  of  the  air/sea  interface,  gives  rise  to  surface  reverberation.  This  type  of 
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reverberation  varies  with  wind  speed  and  transmission  frequency561018.  In  the 
ocean  body,  air  bubbles,  suspended  sediment,  thermal  inhomogeneities,  fish, 
plankton,  and  other  biological  sources  are  the  main  contributors  to  what  is 
classified  as  volume  reverberation11.  Bottom  reverberation  is  caused  by  energy 
scattered  from  the  sea  floor,  and  is  highly  dependent  on  the  floor  type.  Both 
particle  size  and  bottom  relief  are  important  factors37  8.  In  general,  ocean 
acoustic  reverberation  has  a  very  complex  nature,  and  a  highly  variable  power 
spectrum.  In  some  situations,  where  the  sonar  is  fairly  close  to  the  ocean 
boundaries,  or  when  its  beam's  sidelobes  are  pointing  towards  those  boundaries, 
the  energy  reflected  off  of  these  boundaries  makes  a  significant  contribuuon  to 
the  range-DoppIer  map.  This  contribution  usually  has  a  sudden  onset,  may  appear 
at  nonzero  Doppler  frequencies,  and  therefore  may  mask  legitimate  sonar  echoes. 

1.1.  The  statistical  and  spectral  properties  of  reverberation: 

1.1.1.  Amplitude  distribution 

Common  reverberation  models  assume  that  the  reverberation  signal  is 
comprised  of  a  combination  of  numerous  individual  sources  which  scatter  the 
acoustic  energy.  These  sources  scatter  the  energy  with  different  magnitudes  and 
random  phase.  An  expression  for  the  intensity  distribution  of  numerous  sources 
of  equal  magnitude  and  random  phase  has  been  derived  by  Rayleigh51.  The 
probability  that  the  resultant  intensity  will  exceed  the  value  /  is  given  by: 

/ 

P  =  exp  -  —  (1.1) 

/ 

where  I  is  the  average  intensity.  The  actual  received  intensity  is,  of  course,  a 
combination  of  sources  of  many  different  magnitudes,  because  the  individual 
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scatters  are  of  different  size  and  strength,  because  the  projector  may  be 
directional,  and  because  the  transmission  loss  in  different  parts  of  the  ocean  may 
vary.  Nevertheless,  it  can  be  shown  that  Equation  (1.1)  is  still  valid  provided  that 
there  are  a  large  number  of  sources  of  each  magnitude,  and  that  the  number  of 
sources  of  each  magnitude  remains  constant. 

The  applicability  of  the  Rayleigh  model  has  been  tested  in  numerous 
studies44  36  35  ,  for  all  reverberation  classes  (surface,  volume  and  bottom).  All  of 
the  studies  have  found  that  in  general,  the  actual  probability  distribution  of 
reverberation  intensity  fits  the  Rayleigh  model  well.  Moreover,  it  has  been  found 
that  the  received  envelope  of  the  reverberation  is  a  complex  Gaussian  random 
process  (i.e.  it’s  real  and  imaginary  parts  are  jointly  Gaussian  distributed). 

1.1.2.  Coherence 

Coherence  includes  both  the  autocorrelation  and  crosscorrelation  properties 
of  a  statistical  process.  The  observed  reverberation  envelope  of  pulsed  sonar 
possesses  a  "blobby"  nature,  where  the  "biobs"  are  about  a  pulse  width  long.  The 
temporal  autocorrelation  coefficient  of  reverberation  envelope  has  been  found  to 
almost  diminish  in  a  time  delay  equal  to  one  pulse  duration38.  Spatial  correlation 
of  reverberation  received  at  two  vertically  separated  hydrophones  plays  an 
important  role  when  implementing  reverberation  rejection  algorithms.  This 
correlation  has  been  studied,  modeled4  and  measured58'59-57  for  both  volume  and 
boundary  reverberation.  Both  the  models  and  the  experimental  results  show  that 
boundary  reverberation  is  more  correlated  than  volume  reverberation,  and  that 
the  crosscorrelation  coefficient  of  both  types  diminishes  with  increased  frequency 
and  hydrophone  separation. 
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1.1.3.  Frequency  distribution 

The  reverberation  return  resulting  from  a  single  frequency  sonar  pulse 
typically  occupies  a  spectral  band  which  is  non-congruent  to  the  original  pulse 
band.  The  reverberation  band  is  shifted  in  frequency  due  to  Doppler  effects 
caused  by  the  sonar  platform  movement,  combined  with  any  velocities  possessed 
by  the  scattering  media  be  it  volume  scatterers  or  surface  waves.  Furthermore,  a 
spectral  spread  is  observed.  Naturally,  the  finite  duration  of  the  transmission 

1 

pulse,  results  in  a  basic  spread  of  — ,  where  r  is  the  pulse  duration.  Further  spread 

T 

is  a  result  of  the  fact  that  different  Doppler  shifts  are  caused  by  the  reverberation 
returns  received  from  various  directions.  The  Doppler  shift  depends  linearly  on 
'ne  sonar  platform's  velocity  component  in  that  direction,  and  therefore 
reverberation  coming  in  from  different  directions  suffer  different  Doppler  shifts. 
A  third  contributor  to  the  spectral  spread  is  the  one  caused  by  the  random  motion 
of  the  scattering  media. 

A  typical  range-Doppler  map  of  a  reverberation  reti  rn  is  shown  in  Figure 
1.1.  it  was  generated  using  a  reverberation  simulation  package  (REVSIM)  developed 
for  this  research,  which  will  be  described  in  detail  in  a  later  chapter  of  this 
dissertation.  Here,  the  vehicle  was  purposefully  given  a  relatively  high  speed  (30 
m/s),  so  as  to  well  separate  boundary  reverberation  from  volume  reverberation  in 
the  frequency  domain,  and  make  each  of  them  well  defined.  The  range-Doppler 
map  has  also  been  left  shifted  to  compensate  for  the  vehicle's  own  speed.  The 
sonar  array  is  placed  at  a  depth  of  100m,  where  the  ocean  depth  is  400m.  The 
transmit  beam  illuminates  a  sector  ±73’  in  elevation,  and  the  receive  beam  spans 
±60’.  Early  in  the  ping,  around  0.1  sec,  the  onset  of  the  surface  reverberation  is 
well  noticeable.  It  slowly  creeps  towards  0  Doppler  joining  the  volume 
reverberation,  when  aioun  u  0.5  SCC  inC  of  bottom  reverberation  appears.  It 
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too,  then,  slowly  joins  the  volume  reverberation  until  around  1  sec  the  two 
spectral  peaks  become  very  close. 


2.  Element  level  reverberation  time  series  synthesis 


2.1.  Introduction 

The  ocean’s  boundaries  (surface  and  bottom)  and  particulate  matter  in  the 
ocean  volume  scatter  energy  from  an  active  sonar  ping  back  to  the  transducer. 
The  characteristics  of  this  reverberation  are  highly  non-stationary.  A  commonly- 
used  display  of  the  time-evolving  frequency  content  of  reverberation  is  the 
range-Doppler  map.  The  expected  range-Doppler  map  is  useful  in  many  instances 
for  active  sonar  system  performance  characterization. 

As  depicted  in  Figure  2.1.  one  approach  to  estimating  the  expected  range- 
Doppler  map  is  to  propagate  the  pulse  outward  and  calculate  the  expected 
backscattered  return  at  successive  range  increments22.  At  each  range,  the  shell 
representing  the  ensonified  portion  of  the  ocean  and  its  boundaries  is  split  into 
cells  subtending  a  small  azimuth  and  elevation  angle  spread.  The  expected 
spectrum  contributed  by  each  ceil  is  summed  across  all  cells  yielding  the  expected 
spectrum  at  that  range  (i.e.  a  slice  of  the  range-Doppler  map  at  a  given  range). 


Figure  2.1.  Reverberation  mode!  geometry 
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Although  the  expected  range-Doppler  map  is  useful  in  many  instances  for 
active  sonar  system  performance  characterization,  often  the  actual  reverberation 
time  series  is  desired.  Such  reverberation  time  series  are  generated  by  REVGEN 
(REVerberation  GENerator)49’50.  REVGEN  is  a  direct  implementation  of  the  point¬ 
scattering  model  of  oceanic  reverberation  developed  by  Faure,  Ol’shevskii,  and 
Middletonls-45  40’4I-'42.43  jhe  returns  from  a  large  number  of  discrete  scatteres 
distributed  randomly  throughout  the  volume  and  on  the  boundaries  are  summed 
coherently  to  obtain  the  (complex  basebanded)  reverberation  time  series  at  the 
output  of  each  hydrophone  element  or  beam.  The  backscattering  coefficient 
(strength)  for  each  reverberation  type  (surface,  volume,  and  bottom)  is  specified 
along  with  random  scatterer  motion,  platform  trajectory,  absorption,  boundary 
reflection  losses,  and  transmitting/receiving  beam  patterns.  REVGEN  has  proved 
very  effective  in  the  simulation  of  oceanic  backscatter  in  a  wide  range  of 
applications7’63'  19-2’3. 

Although  very  good,  REVGEN  also  is  very  computationally  intensive  and  is 
not  suitable  for  generating  a  large  number  of  pings  for  a  given  environment  as 
needed  when  performing  Monte  Carlo  simulations  for  „:tive  sonar  system 
performance  characterization.  Typically,  generating  a  single  sonar  ping  similar  to 
the  one  synthesized  in  Section  2.3,  will  take  REVGEN  a  few  hours  to  run.  In  order 
to  generate  a  large  ensemble  of  pings,  REVGEN  would  have  to  be  run  hundreds 
and  maybe  thousands  of  hours.  An  alternative  to  the  point  scattering  approach 
taken  in  REVGEN  is  to  synthesize  reverberation  by  passing  a  vector  white  noise 
process  through  a  multichannel  HR  filter  whose  time-varying  transfer  function 
matches  the  time-evolving  spectral  characteristics  of  the  expected  auto  and  cross 
range-Doppler  maps.  Of  specific  interest  is  the  generation  of  array  element-level 
or  beam-level  reverberation  time  series  which  have  appropriate  element-to- 
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element  or  beam-to-beam  crosscorrelations.  Related  approaches  are  discussed 
jn20,62.9.34  Obviously,  one  run  time  investment  here  is  in  generating  the 
multidimensional  time-varying  HR  filter  coefficients.  Surprisingly,  this  task  in 
itself  is  relatively  short,  and  typically  takes  a  few  minutes.  Once  the  filter 
coefficients  are  available,  generating  each  sonar  ping  takes  a  few  minutes,  as 
contrasted  to  a  few  hours  using  REVGEN.  The  time  savings  factor  exceeds  100, 
and  allows  the  synthesis  of  a  complete  sonar  ping  ensemble  over  the  course  of  a 
few  hours.  This  makes  the  use  of  a  Monte  Carlo  simulation  approach  feasible. 

Section  2.2  will  describe  one  approach  to  the  generation  of  the  expected 
auto  and  cross  range-Doppler  maps  and  how  to  fit  a  time-varying,  multichannel, 
autoregressive  (AR)  time  series  model  to  these  range-Doppler  maps.  Then,  a  set 
of  simulation  results  is  provided  in  Section  2.3  and  is  compared  with  REVGEN's 
output. 

2.2.  Multichannel  revereration  time  series  synthesis 

2.2.1.  Expected  range-Doppler  map  generation  (RVMD) 

There  are  several  approaches  which  can  be  used  to  calculate  the  expected 
auto  and  cross  range-Doppler  maps.  The  specific  environmental  model  used  as 
the  basis  for  the  simulations  presented  in  Section  2.3  is  a  multibeam  extension  of 
the  approach  discussed  in22.  Array  elements  can  be  viewed  simply  as  beams  with 
broad  spatial  response  characteristics,  when  an  element  level  Range-Doppler  map 
is  sought.  The  module  whose  function  is  to  compute  the  scattering  function 
resulting  from  the  combined  effects  of  the  environment,  vehicle  dynamics  and 
transmit  and  receive  beam  patterns  is  modified  so  it  handles  a  pair  of  receive 
beams  at  a  time.  When  the  program  is  handling  an  auto-scattering  function,  both 
receive  beams  are  identical,  but  when  the  program  is  dealing  with  a  cross- 
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scattering  function,  the  two  receive  beams  are  the  ones  whose  cross  scattering 
function  is  being  computed.  The  programs  run  until  all  possible  receive  beam 
pairs  are  exhausted.  The  result  is  a  multidimensional  scattering  function. 

As  shown  in  Figure  2.1,  the  reverberation  model  geometry  consists  of  a 
spherical  shell  representing  the  portion  of  the  ocean  illuminated  by  the  signal 
wavefront  at  the  instant  in  time  corresponding  to  a  range  R  after  transmit.  The 
expected  range-Doppler  map  is  calculated  in  two  stages. 

First,  the  multidimensional  scattering  function  is  calculated  following  the 
steps  outlined  in  Figure  2.2.  The  illuminated  volume  of  the  ocean  is  divided  into 
cells.  An  evaluation  of  scatterer  motion  relative  to  the  platform  is  made  for  a 
measure  of  spectral  shifting  and/or  spreading  due  to  the  environment.  For  each 
cell,  and  each  beam  pair,  the  sonar  equation  is  solved  in  order  to  determine  the 
gross  attenuation  experienced  by  a  transmitted  signal  during  the  scattering 
process.  The  total  scattering  function  for  each  beam  pair  for  surface,  volume,  and 
bottom  reverberation  is  computed  by  summing  the  incremental  contributions  from 
all  cells.  These  scattering  functions  include  the  effects  of  platform  motion, 
!  transmit  beampattern,  both  receive  beam  patterns,  and  the  environment  (surface, 

j  volume,  and  bottom  backscattering  strengths,  scatterer  velocity  distributions, 

surface  waves  and  current  layers,  and  sound  absorption).  An  isospeea  sound 
speed  profile  is  assumed  and  reflections  at  the  surface  and  bottom  boundaries  are 
not  permitted. 
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Figure  2.2.  Scattering  function  calculation 


In  Figure  2.2,  X  is  a  time  delay  (range)  variable,  is  the  Doppler  shift 

accounting  for  vehicle  motion  relative  to  the  [i,j)th  grid  cell  and  mean  scatterer 

2  V 

speed  relative  to  a  stationary  vehicle,  —  fr  is  the  frequency  offset  such  that  a 

C 

stationary  target  on  vehicle  axis  has  zero  Doppler  shift,  ,  is  the  zero-mean 
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Gaussian  speed  distribution  of  scatterers  in  (i.j)  grid  cell,  5^  ( f,\)  ,  is  the 
Doppler  shifted  speed  distribution  of  scatterers  in  (i,j)th  grid  cell,  Ajjkl  is  the 
scattering  level  (normalized  to  1  s  in  range)  from  the  (i,j)th  grid  cell,  (k,l)th  beam 
pair,  is  the  scattering  function  for  the  (i,j)th  grid  cell,  {k,l)th  beam  pair, 

~k  1 

and  (f,\)  are  the  surface,  volume,  bottom  and  total  scattering  functions  for  the 
(k,l)th  beam  pair. 

Outlined  in  Figure  2.3,  the  second  stage  completes  the  reverberation 
model  by  combining  a  detailed  description  of  the  transmit  signal  (i.e.  pulse  length, 
envelope  shape,  and  source  level)  with  the  characterization  of  the  environment 
provided  by  the  multidimensional  scattering  function.  The  convolution  of  the 
scattering  function  with  the  transmitted  spectrum  yields  the  expected  multi-beam 
range-Doppler  map. 
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Figure  2.3.  Convolution  with  the  transmitted  spectrum 
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In  Figure  2.3  f(t)  is  the  complex  envelope  of  transmit  pulse,  is  the 

normalized  energy  spectrum  of  the  transmit  pulse,  Et  is  the  transmit  pulse  energy, 
Sr’  ( /'.X)  are  the  surface,  volume,  bottom  and  total  scattering  functions  for  the 

th  ~k  l 

( k,l )  beam  pair,  and  are  the  surface,  volume,  bottom  and  total 

reverberation  power  spectra  for  the  ( k,l)th  beam  pair. 

2.2.2.  Time-varying  IIR  filter  derivation 

The  classical  linear  minimum  mean  square  error  (MMSE)  estimation 
problem  has  been  developed  by  Kolmogorov29  and  Wiener68,  and  has  been  since 
one  of  the  foundations  of  estimation  theory.  It  is  discussed  is  standard  estimation 
theory  texts34.  Its  development  is  therefore  left  in  this  dissertation  for  an 
Appendix  (A). 

The  end  result  of  the  linear  MMSE  problem  are  the  well  known  normal 
equations,  in  this  case,  their  multichannel  version: 


AqRq+AoJH- i+ . . . +APR~P  —  G0  (2.1) 

AqK\+AqRq+ . +APR_P+ 1  =  G  i 

AqRp  +Ap_  i  o+ . +ApRq  =  Gp 


or  in  their  matrix  form: 


[40^! . 


Ro  Ri 
R-\  Rq 


p  1 


Rp 

Rp-x 


G0£{ . xGp 


(2.2) 


Ro 


Since  we  are  dealing  with  an  N-input  N-output  system.  A,  are  all  NxN  matrices,  and 
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so  are  Rj  and  Also,  are  the  crosscorrelation  matrices  between  the  input  and 
the  desired  signal.  In  the  specific  problem  of  linear  prediction,  those  reduce 
simply  to  the  crosscorrelation  between  the  input  at  time  n  and  the  input  at  time 
(n- 1),  or  the  input  autocorrelation  at  lag  (f+1): 

£,  =  Ri+ 1  (2.3) 

In  this  case,  the  matrix  form  of  the  normal  equations  takes  the  form: 


R0  Rx 


R 


-l 


Rp 

Vi 


(2.4) 


R-.  R-> 


■P+ 1 


R0 


Solving  the  above  equations  will  yield  the  multichannel  filter  coefficients 

Our  goal  is  to  pass  a  vector  white  uncorrelated  noise  process  through  a 
multichannel  HR  filter  as  shown  in  Figure  2.4  yielding  a  vector  time  series  at  the 
output  which  approximates  the  time-evolving  (auto  and  cross)  spectral 
characteristics  of  reverberation  as  it  would  be  seen  at  the  output  of  the  elements 
or  beams  of  a  transducer  array.  The  time-varying  coefficients  of  the  multichannel 
one  step  forward  prediction  filter  can  be  derived  in  a  number  of  ways. 
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a(n) 


Figure  2.4.  Autoregressive  process  generation  model 


Here,  slices  at  a  given  range  from  the  set  of  expected  auto  and  cross 
range-Doppler  maps  will  be  viewed  as  the  multichannel  power  spectrum  of  a 
stationary  vector  random  process.  First,  the  multichannel  power  spectrum  is 
inverse  Fourier  transformed  to  yield  the  corresponding  multichannel  correlation 
function.  Then,  a  multichannel  extension  of  the  Levinson-Durbin  algorithm  is 
used  to  derive  the  one-step  forward  prediction  error  (inverse)  filter  shown  in 
Figure  2.5  which  provides  the  HR  filter  coefficients.  The  Levinson-Durbin 
algorithm  has  been  discussed  extensively  in  numerous  publications33-14-65  69  52  28, 
and  its  development  will  not  be  repeated  here.  It  takes  advantage  of  the  special 
structure  of  the  autocorrelation  matrices,  to  achieve  an  efficient  method  of  matrix 
inversion. 
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Figure  2.5.  One  step  forward  prediction  error  filter 


Each  increment  in  range  is  processed  in  exactly  the  same  way  yielding 
successive  sets  of  HR  filter  coefficients.  The  prediction  error  filter  order  is 
chosen  such  that  it  is  of  the  smallest  order  large  enough  to  accomplish  emulating 
ah’  of  the  spectral  contents  in  the  input  data.  In  the  type  of  problems  encountered 
in  this  dissertation,  where  the  input  spectrum  typically  posesses  two  spectral 
humps,  namely  the  main  volume  reverberation  spectral  hump,  and  anothe. 
boundary  reverberation  spectral  hump,  experience  shows  that  filter  order  can  be 
kept  to  a  low  value  (3-4),  while  preserving  all  spectral  features  of  the  input.  This 
is  well  demonstrated  in  the  simulations  in  Section  2.3. 

2.2.3.  Reverberation  time  series  synthesis 

The  sets  of  HR  filter  coefficients  (one  set  per  range  increment)  then  are 
used  as  the  coefficients  of  a  time-varying,  multichannel  HR  filter.  Since  each 
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coefficient  set  represents  a  range  slice,  the  range  slices  originally  selected  for 
RVMD  have  to  be  narrow  enough  so  that  the  time  varying  spectra  are  relatively 
smooth  from  range  slice  to  the  next.  A  vector  white  uncorrelated  noise  process  is 
passed  through  the  filter  yielding  a  vector  output  time  series  which  simulates  the 
outputs  of  a  real  transducer  array  in  an  ocean  whose  environmental  parameters 
are  those  input  to  RVMD. 

2.3.  Active  sonar  simulation 

As  an  example  of  this  approach  to  reverberation  time  series  synthesis, 
consider  the  following  active  sonar  system.  An  array  of  9  rows  is  mounted  on  a 
submerged  vehicle  as  shown  in  Figure  2.6. 

1  surface 


_ y  bottom 

Figure  2.6.  Submerged  Array  Example 

The  vehicle  is  traveling  at  25  knots  at  a  depth  of  30  m.  Ocean  depth  is  180  m.  The 
sonar  transmits  a  180  ms,  rectangularly  windowed  pulse. 
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The  9  row  array  has  a  sensor  spacing  of  X/2  at  the  operating  frequency  of 
the  sonar.  A  conventional  beam  is  formed  by  simply  summing  the  outputs  of  all 
array  rows.  The  corresponding  beam  pattern  and  a  polar  plot  of  a  vertical  slice 
through  the  beam  pattern  at  0“  bearing  are  shown  in  Figure  2.7. 

A  difference  beam  is  formed  by  subtracting  the  outputs  of  the  2  center  rows  of 
the  array.  The  corresponding  beam  patterns  similar  to  those  in  Figure  2.7  are 
shown  in  Figure  2.8. 
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A  transmit  beam  is  assumed  which  uniformly  illuminates  the  medium  over  a  sector 
±30°  in  bearing  and  ±75°  in  elevation  normal  to  the  array  plane. 

The  backscattered  return  to  this  sonar  was  simulated  using  REVGEN  and 
the  multichannel  HR  filter  approach  discussed  in  Section  2.2.  In  this  case,  an 
isovelocity  sound  speed  profile  was  assumed  (c  -  1500  m/s)  along  with  the 
following  parameters:  (a)  volume  backscattering  coefficient,  sy  -  -70  dB,  (b) 
bottom  backscattering  coefficient,  s^  -  -20  dB,  and  (c)  surface  backscattering 
coefficient,  s$  =»  -30  dB,  and  (d)  attenuation  due  to  absorption,  a  -  4  dB/km. 

REVGEN  range-Doppler  maps  of  the  returning  reverberation  for  the 
conventional  (sum)  and  difference  beams  as  well  as  the  the  corresponding  cross¬ 
power  range-Doppler  map  are  shown  in  Figure  2.9  (fs  -  1.27  kHz).  These  were 
generated  by  taking  successive  128-point  FFT’s  (Kaiser-Bessel  window,  a-2.5) 
overlapped  by  87.5%  (16  points).  In  addition,  the  range-Doppler  maps  have  been 
left-shifted  to  compensate  for  platform  velocity.  In  Figure  2.9a  (sum  beam), 
volume  reverberation  dominates  early  in  the  range-Doppler  map.  Later  (at 
approximately  0.25  s  in  range),  the  onset  of  bottom  reverberation  arriving  through 
the  side  lobes  of  the  main  beam  is  seen,  in  Figure  2.9b  (difference  beam),  volume 
reverberation  also  dominates  early  in  the  range-Doppler  map.  The  onset  of  bottom 
reverberation  is  more  pronounced  here  than  in  Figure  2.9a  due  to  the  response 
characteristics  of  the  difference  beam  pattern. 
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Figure  2.9. 


REVGEN  -  Sum,  Difference  and  Cross  range-Doppler  maps. 


The  multichannel  HR  filter  results  will  be  presented  following  the  three 
steps  discussed  in  Section  2.2.  The  RVMD  parameters  were  identical  to  those  used 
in  the  REVGEN  simulation.  First,  Figure  2.10  shows  the  RVMD  auto  and  cross 
range-Doppler  maps.  Note  the  "scallops"  in  the  surface  portion  of  the  map,  which 
are  due  to  the  sidelobes  intersecting  the  ocean  boundary.  Second,  the 
corresponding  time-evolving  inverse  linear  predictor  spectral  estimates  (filter 
order  p-4)  are  shown  in  Figure  2.11.  They  indicate  how  closely  the  derived  HR 
filter  matches  the  expected  auto  and  cross  range-Doppler  maps.  Lastly,  Figure 
2.12  shows  the  results  of  passing  a  vector  white  noise  process  through  the 
multichannel  HR  filter  and  processing  the  output  as  done  with  the  REVGEN  time 
series.  As  is  easliy  seen,  the  results  are  quite  similar  to  the  REVGEN  range-Doppler 
maps  in  Figure  2.9. 

Figure  2.13  presents  expanded  high-resolution  range-Doppler  maps  of  the 
main  beam,  in  order  to  further  emphasize  the  similarity  of  REVGEN  output  (2.13a), 
RVMD  output  (2.13b),  the  linear  predictor  (2.13c),  and  finally  REVSIM  output 
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Figure  2.1 1.  Inverse  Linear  Predictor  -  Auto  and  Cross  range-Doppier  maps. 
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Figure  2.12.  REVSIM  -  Auto  and  Cross  range-Doppler  maps 
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Figure  2.13.  Main  beam  •  high  resolution  range-Doppier  maps 


3.  Signal  detection  in  boundary  reverberation 


As  described  in  a  previous  chapter,  ocean  acoustic  reverberation  has  a 
very  complex  nature,  and  a  highly  variable  power  spectrum.  Reverberation  is  a 
source  of  interference  which  active  sonar  systems  have  to  combat.  In  some 
situations,  where  the  sonar  array  is  fairly  nose  to  the  ocean  boundaries,  the 
energy  reflected  off  of  these  boundaries  (bottom  cr  surface),  makes  a  significant 
contribution  to  the  range-Doppler  map.  This  contribution  usually  has  a  sudden 
onset,  may  appear  at  nonzero  Doppler  frequencies,  and  therefore  may  mask 
legitimate  sonar  echoes.  The  shape  of  the  transmitting  and  receiving  beams  of  the 
active  sonar  modify  the  time  evolving  range-Doppler  map  observed  by  the 
receiver.  Traditionally,  sonar  systems  have  discriminated  against  boundary 
reverberation  by  forming  fixed  receiving  beams  which  had  low  sidelobe 
characteristics  in  the  direction  of  the  ocean  boundaries.  However,  the  energy 
leaking  through  the  low  sidelobes  may  still  be  a  major  contributor  to  the  overall 
noise  background  level.  Furthermore,  since  the  sonar  system  is  typically  moving 
and  varying  its  depth  in  the  water  column,  the  direction  of  boundary 
reverberation  is  not  fixed. 

The  above  reasoning  led  many  researchers  to  propose  using  an  adaptive 
structure  in  order  to  track  the  interference  direction  and  to  place  spatial 
beamformer  nulls  in  that  direction.  Such  a  beamformer  may  continuously  learn  the 
boundary  interference  direction  and  adjust  itself  to  cancel  it  out,  thereby 
enhancing  signal  to  noise  ratio  and  improving  the  detection  performance. 
Implementation  of  such  an  ad-hoc  structure  follows  the  lines  of  the  well  known 
adaptive  noise  canceler66.  In  addition  to  a  main  beam  which  receives  well  in  the 
desired  look  direction,  one  or  more  reference  beams  are  formed.  The  reference 
beams  receive  well  in  the  direction  of  the  boundary  interference,  and  have  spatial 
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nulls  in  the  main  look  direction.  The  output  of  the  main  beam  contains  both  the 
desired  signal,  and  a  contaminant  which  is  the  contribution  of  the  boundary 
interference  leaking  through  side  lobes  or  the  edges  of  the  main  lobe.  Ideally,  the 
reference  beams  contain  only  a  replica  of  the  interference.  The  output  of  the 
reference  beam,  or  beams,  is  processed  by  an  adaptive  filter  and  then  subtracted 
from  the  main  beam.  The  adaptive  filter  tries  to  provide  a  good  estimate  of  the 
interference  portion  of  the  main  beam  output,  and  the  final  error  output  ideally 
contains  only  the  desired  signal.  Implementation  of  such  adaptive  reverberation 
cancellation  schemes  is  reported  by  Hodgkiss  and  Alexandrou23.  The  boundary 
reverberation  canceled  output  of  the  adaptive  structure  can  now  be  treated  as 
containing  a  known  signal  in  noise  (the  volume  reverberation),  a  classical,  solved 
problem. 

The  proposed  solution  discussed  so  far  is  representative  of  an  approach 
which  is  based  on  intuition.  A  typical  conventional  adaptive  processor  would 
combine  one  or  more  of  the  following  building  blocks  in  order  to  make  a  decision 
about  target  presence: 

1.  Fixed  Beamformer. 

2.  Adaptive  beamformer. 

3.  Adaptive  nullformer. 

4.  Adaptive  spectral  whitener. 

5.  Matched  filter. 

6.  Threshold  detector. 

Use  of  each  of  the  above  building  blocks  is  intuitively  reasonable,  but  it  is  not  at 
all  clear  that  any  combination  of  these  :omponents  yields  an  optimal  processor 
(under  any  optimization  criterion). 

Another  approach  suggests  treating  the  problem  as  a  whole  right  from  the 
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beginning  without  imposing  intuitive  components  on  the  processor  structure,  and 
using  all  a-priori  knowledge  available.  Detection  theory  provides  us  with  a 
mathematical  framework  out  of  which  optimum  processors  can  be  designed54'64. 
The  processor  will  evolve  out  of  the  mathematical  solution  of  the  problem,  and 
will  not  be  restricted  to  using  familiar  structures. 

Although  Bayes  optimal  processors  have  been  derived  for  the  case  of 
volume  reverberation55,  little  work  has  been  done  which  takes  advantage  of  a- 
priori  knowledge  of  the  time-evolving  spatial  characteristics  of  boundary 
reverberation.  Related  Bayes  optimal  work  concerning  interference  sources  of 
certain  and  uncertain  (but  not  time  varying)  location  is  contained  in1  and  26-24'25 
respectively. 

Here,  a  classical  detection  theoretic  approach  is  applied  to  the  processing 
of  a  vector  time  series.  That  vector  may  be  composed  of  the  single  array  element 
outputs,  or  of  some  preformed  beams.  An  optimality  criterion  is  chosen,  and  then 
the  processor  structure  is  allowed  to  evolve  freely  out  of  the  mathematical 
solution  of  the  problem.  Any  uncertain  parameters  are  treated  as  random  variables 
and  all  knowledge  about  them  is  summarized  in  a-priori  probability  density 
functions. 

3.1.  The  adaptive  (noise  canceling)  ad-hoc  detector 

As  previously  mentioned,  the  time  evolving  nature  of  ocean  acoustic 
boundary  reverberation  has  led  many  researchers  to  try  applying  adaptive 
filtering  schemes  to  an  observed  vector  time  series,  in  order  to  cancel  out,  or 
reduce  the  interference  level  present  in  it.  The  idea  of  adaptively  canceling 
interference  sources  thereby  enhancing  the  desired  signal  is  intuitively 
reasonable,  and  seems  very  attractive.  It  is  not  at  all  clear,  though,  that  using  such 
a  canceler  as  a  building  block,  following  it  by  another  ad-hoc  detection  building 
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block  is  a  globally  optimal  solution  from  any  optimal  detection  criterion  point  of 
view. 

As  is  implied  from  its  name,  an  adaptive  filter  "adapts"  to  the  changing 
conditions  by  varying  its  weights,  and  going  through  a  learning  process.  Various 
mechanisms  that  control  the  variation  of  the  filter  weights  have  been  proposed 
and  implemented.  Each  one  of  these  mechanisms  is  based  on  some  optimality 
criterion  (e.g.  minimum  mean  square  error),  which  is  local  to  the  adaptive  filter 
structure,  and  does  not  necessarily  conform  to  any  global  optimum  detection 
criterion. 

The  adaptive  filtering  structures  may  be  viewed  as  realizable  approximate 
solutions  to  the  classical  optimal  Wiener  filtering  problem,  and  can  be  developed 
both  from  a  statistical  point  of  view  (assuming  stationarity  and  ergodicity),  or  by 
using  a  deterministic  least  squares  approach.  Both  approaches  can  be 
implemented  through  block  processing  or  time-recursive  techniques,  and  may  be 
implemented  in  direct  form,  or  through  efficient  lattice  structures.  Built  into  all  of 
the  adaptive  filter  implementations  there  is  an  adaptation  coefficient  which  serves 
as  a  ’forgetting  factor’  (i.e.  weighting  less  past  in  favor  of  more  recent  data). 

There  are  three  general  classes  of  adaptive  structures.  The  first  class  is 
based  on  the  method  of  steepest  descent  and  is  generally  referred  to  as  the  LMS 
algorithm.  The  work  of  Widrow  and  Hoff  67  has  served  as  the  basis  to  this  method. 
The  second  class  is  the  gradient  lattice  (GL),  and  the  third  is  of  the  deterministic 
least  squares  type,  and  is  referred  to  as  the  least  squares  lattice  (LSL).  A 
comparison  between  these  three  classes  is  beyond  the  scope  of  this  work,  and 
Alexandrou4  is  a  good  reference  in  that  regard.  Due  to  its  superior  behavior 
during  abrupt  changes  in  the  processed  data,  the  LSL  class  of  adaptive  filters  has 
been  typically  selected  by  researchers  dealing  with  the  ocean  boundary 
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reverberation  problem.  Hodgkiss  and  Alexandrou23  and  Alexandrou3  have  shown 
how  the  LSL  filter  can  be  used  to  cancel  or  reduce  sea  surface  reverberation 
interference.  LSL  class  adaptive  filters  are  used  in  this  dissertation  as  part  of  an 
ad-hoc  detection  scheme.  In  fact,  a  new  structure  of  a  joint  process  pole-zero 
adaptive  LSL  filter  is  developed  in  the  next  chapter,  and  its  performance  is 
compared  to  that  of  the  standard  all-zero  LSL  filter. 

The  ad-hoc  detection  scheme,  then,  accepts  as  inputs  a  main  beam  output 
which  contains  the  signal  contaminated  by  interference,  and  one  or  more 
reference  outputs,  which  contain  the  interference  alone.  The  adaptive  filter  tries 
to  cancel  the  interference,  and  the  filter's  output  is  then  presumed  to  be 
interference  free,  and  is  presented  to  a  signal  known  exactly  in  noise  known 
exactly  Bayes  optimum  detector  for  making  the  decision.  Figure  3.1  illustrates  the 
structure  of  the  ad-hoc  detector. 


Figure  3.1  The  ad-hoc  detector 
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3.2.  Complex  exponential  Fourier  series  representation  of  signals 

In  order  to  be  able  to  use  classical  statistical  techniques  and  apply  those  to 
optimal  detection  theory,  we  need  to  obtain  a  finite  dimensional  observation 
vector  which  is  then  used  in  computing  the  joint  probability  density  functions 
under  the  true  and  the  null  hypotheses.  These  probability  density  functions,  in 
turn,  are  used  to  form  the  likelihood  ratio.  Typically,  the  sector  of  observables 
available  for  processing  is  that  of  the  array  element  outputs,  or  some  preformed 
beam  outputs.  In  any  case,  this  is  a  vector  of  continuous  time  waveforms,  which 
means  that  some  mapping  technique  is  required  in  order  to  collapse  this  space 
into  a  finite  dimensional  space.  There  are  three  well  known  techniques  to  perform 
this  transformation,  namely  time  sampling,  Karhunen-Loeve  expansion,  and 
trigonometric  Fourier  series.  The  best  selection  of  the  mapping  technique  usually 
depends  on  the  problem  one  is  faced  with,  and  a  good  choice  can  often  reduce  the 
complexity  of  writing  the  likelihood  ratio  expressions. 

Adams  and  Nolte1  have  shown  that  using  the  trigonometric  Fourier  series 
mapping  technique  leads  to  significant  mathematical  tractability,  and  Hodgkiss27 
derived  the  conditions  that  have  to  be  met  in  order  to  be  able  to  represent  the 
observed  vector  in  terms  of  its  Fourier  coefficients.  Certainly,  with  the  advent  of 
fast  computers,  and  «.ne  availability  of  fast  Fourier  transformation  algorithms, 
using  this  approach  seems  very  attractive.  The  mathematical  derivations  of 
optimal  Bayes  detectors  in  this  dissertation  follow  aimilcu'  derivations  made  by 
Adams  and  Nolte1  and  Hodgkiss27  except  that  their  derivations  dealt  with  the 
passive  sonar  problem,  where  this  work  is  dealing  with  its  active  counterpart.  The 
trigonometric  Fourier  series  mapping  is  used,  then,  throughout  this  work. 

T  T_ 

2  ’  2  , 


It  is  well  known  that  a  well  behaved  function  f(t)  on 


can  be 
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represented  as: 


f(t)  =  lim  £>(n)  —  exp(jnu0t) 


N—  ao_N  l  7 


where  =  — ,  and 


/  s1^  Tf2 

1 

a(n)  =  —  I  f(t)exp(-jnuQt)dt 
l'T'  -772 


If  fit)  is  bandlimited  and  has  no  dc  component  then  on  f(t)  can  be 

12  2 


expressed  as: 


fit)  =  2  Re  £>(n)  —  exp ijnu0t) 

w,  T> 


Then  f{t)  is  mapped  into  a  finite-dimensional  vector/. 


/T  »  l«(JV1).«(Ni+ 1) . fl(Nz)]  (3.4) 

The  above  can  also  be  shown  for  a  stationary  random  process  rit),  with  power 

T  7 

spectral  density  R(w).  It  can  be  expressed  on  - —  , —  as: 

2  2 


1 

rit)  =  lim  J]z(n)(— )1/2expO'nw0f) 
N— ►  00  ..  7 


where  z(n)  is  a  random  variable  defined  as  in  (3.2)  above. 


If  r(t)  is  bandlimited  and  has  no  dc  component  then: 


rit)  =  2  Re  £z(n)(— )l/2exp(jnu,0t) 

n,  T 
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i.e.  every  realization  of  the  random  process  r(t)  can  be  mapped  into  a  vector: 

zT  =  [z(N1),z(W1+l) . z(n2)  ]  (3.7) 

If  r(t)  is  a  Gaussian  random  process,  then  z(n),  which  is  a  linear  functional 
of  k{()  is  a  Gaussian  random  variable.  For  zero  mean  noise  processes,  z  has  the 
zero  vector  as  its  mean.  For  large  7,  Papoulis47  has  shown  that  the  components  of 
(3.7)  are  approximately  uncorrelated  and  z(n)  has  a  variance  R(nw0).  The 
covariance  matrix  of  z  is  therefore  diagonal.  It  can  also  be  shown  that  2  is  a 
complex  Gaussian  vector,  i.e.  is  a  (N2  -  N,  +  l)-tuple  of  complex  random  variables 
such  that  the  vector  of  real  and  imaginary  parts  is  (N2  -  N,  +  l)-variate  Gaussian. 

3.3.  The  array  detection  problem 

An  array  of  receiver  elements  observes  a  vector  of  time  waveforms  as 

r  1 1  _ _ 

«.  GtiO  w  o. 

rr(t)  =  lro(07,(t) . rK_,(t)]  (3.8) 

where  subscript  denotes  the  array  element. 

If  all  the  received  waveforms  are  bandlimited  then  according  to  (3.2)  one  can 
express  them  in  a  Fourier  expansion,  as  follows: 

'f 

-  exp  (jnw0r)  (3.9) 

The  Fourier  coefficients  are: 

1/2  772 

z.k{n)  =  —  /  rk(t)exp(-jn±j0t)dt  (3.10) 

T’  -t/ 2 

One  can  now  group  the  Fourier  coefficients  for  a  single  frequency  index  but  from 


N 

rk{t)  =  lim  £zk(n) 

N—  00  ,, 
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all  the  array  elements  in  a  K  dimensional  vector  as  follows: 

ZT{n)  =  (z0(tt).z,(n) . zK  _  [(«)]  (3.1 1) 

and  group  now  vectors  from  all  the  frequency  indices: 

ZT  =  lZT(Nl)£T(Nl  +  1) . ZT(N2)]  (3.12) 

In  this  way  the  time  waveform  observed  on  the  K  elements  is  mapped  into  a 
K{N2  -  N !  +  1)  dimensional  vector.  It  can  be  shown  that  Z  is  a  complex  Gaussian 
vector,  and  that  its  covariance  matrix  (for  large  T),  is  block  diagonal.  There  are 
(N2  -  N j  +  1)  KxK  blocks,  each  representing  the  correlation  between  array 
elements  at  a  specific  frequency. 

3.4.  Hypothesis  testing 

The  theory  of  hypothesis  testing  began  with  the  work  of  Bayes6  on 
conditional  probability.  The  notions  of  cost  and  risk  were  introduced  by  Wald62, 
who  initiated  much  of  the  recent  theoretical  work  along  the  lines  of  the  theory  of 
games.  Communications  receivers  based  on  the  ideas  of  conditional  probability 
were  proposed  by  Kotel’nikov30  in  a  dissertation  from  1947,  published  in  1956 
and  translated  into  English  in  1959.  At  the  same  time  Woodward  and  Davies70, 
and  Woodward71  suggested  applying  conditional  probability  to  signal  detection. 
An  analysis  of  signal  detection  within  the  framework  of  conditional  probability 
was  performed  by  Middleton39,  and  Peterson,  Birdsall  and  Fox  48  were  at  the  same 
time  designing  receivers  based  on  the  likelihood  ratio.  The  following  summarizes 
the  thought  process  developed  by  the  above  mentioned  researchers. 

A  processor  which  has  available  an  observed  vector  time  series  z.  has  to  be 
designed.  The  processor  is  forced  to  make  a  binary  decision  between  two  mutually 
exclusive  and  exhaustive  hypotheses,  namely,  signal  absent  or  signal  present: 
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H0:z  =  ii(20)  (3.13) 

Hj  :  z  =  .$(£,)  +  ndO  (3.14) 

where  a(^),  s(S.{)  and  ndi)  are  noise  and  signal  vectors  written  as  functions  of 
paiameter  vectors  Sc  and  lv  The  parameter  vectors  determine  the  shape  and 
character  of  the  noise  and  signal  vectors,  and  may  be  fully  known,  partially  known, 
or  even  completely  unknown. 

The  processor’s  choice  results  in  one  of  four  different  outcomes,  i.e  a  false 
alarm,  a  miss,  a  detection  or  a  null  decision.  Their  corresponding  probabilities  are: 


False  Alarm:  Q10  =  / p0(z)dz  (3.15) 

R  1 

Miss:  Qo,  =  /  p{(z)dz  (3.16) 

RO 

Detection:  Qu  =  /  Pi(z)dz  (3.1 7) 

R1 

Null  Decision:  Qoo  =  /  p0(z)dz  (3.18) 

RO 


where  p0(z)  and  p^z)  are  the  probability  density  functions  of  the  observed  vector 
under  the  two  hypotheses. 

A  cost  is  attached  to  every  one  of  the  outcomes,  and  consequently  a  cost 
matrix  is  defined: 


C  = 


(3.19) 


Where  CtJ  is  the  cost  of  choosing  when  Hj  is  true.  The  cost  selection  depends 
heavily  on  the  scenario  of  the  problem,  but  usually  the  costs  for  making  a  correct 
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decision  are  chosen  lower  than  those  for  an  incorrect  decision.  A  negative  cost 
means  »  reward  (conceivable  for  a  correct  decision). 

A  risk  function  is  also  attached  to  every  outcome,  and  is  defined  as  follows: 

R-CjjQjj  (3.20) 

This  function  is  the  product  of  the  cost  and  the  probability  of  an  outcome,  and 
represents  the  risk  incurred  in  case  that  outcome  materializes. 

3.5.  The  Bayes  Criterion 

ihe  tsayes  criterion  attempts  to  minimize  the  average  risk,  and  seems  a 
natural  approach  for  an  observer  forced  to  make  a  large  number  of  decisions 
under  repetitive  circumstances.  If  $■  is  the  probability  that  H0  is  true,  the  average 
risk  is: 


R  —  flQjoQoo  +  ^ioQioI  +  (1  —  f)[C01Qoj  +  CUQU]  (3.21) 

=  dC00J P0{z)dz  +  Ci0Jp0(z)dz]  +  (1  -  ?)( C0ifpi(z)dz  +  CnJp,(z)dz] 

R0  R1  Ro  R1 

Now,  since  / Pj(z)dz  =  1  -  / Pt{z)dz  /=0,1 

*o  R. 

Then 


R  =  flCoo  +  (CI0  —  Cqo  )/ Po(z)dz\  +  (1  —  f)[  Cq  i  +  (C0,  —  Cu)Jp1(z)dz]  (3.22) 


=  fQo  +  (1  -  f)C01  -  (1  -  f)(C0,  -  Cn)f[p,(z)  -  A 0p0(z)]dz 


Where: 
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?((- io  ~  C00) 

(1  -f)(C01  -  Cu) 


(3.23) 


Furthermore,  since  (C01  -  Cn)  is  positive,  minimizing  R  means  to  select  Rx  such 
that  the  integral  is  as  large  as  possible.  In  other  words,  choose  K,  such  that  all 
points  for  which  the  expression  px(z)  -  A0p0Cz)  is  positive  are  included. 


This  leads  to  the  following  decision  rule: 


P  iCz) 

=  Lz  :  ACz)  = -  >  A0j  (3.24) 

PoU) 

P\U) 

R0  =  [z  .  A (z)  =  -  <  A0] 

PoU) 

The  optimum  Bayes  detector  in  the  least  risk  sense  is  implemented,  therefore,  in 

PXU) 

forming  the  ratio  A (z)  = -  which  is  called  the  likelihood  ratio.  The 

PoU) 

likelihood  ratio  is  compared  to  a  prescribed  threshold,  and  the  decision  is  made 
based  on  the  comparison's  outcome.  The  above  procedure  serves  as  the  basis  for 
all  the  optimum  detector  derivations  in  this  dissertation. 

3.6.  Implementing  the  Bayes  criterion  for  the  detection  of  a  known  signal 
in  spatially  correlated  Gaussian  noise 

The  following  outlines  the  procedure  used  when  deriving  Bayes  optimum 
detectors.  This  case  is  a  very  simple  one,  in  which  everything  is  known  about  the 
signal  and  the  noise,  and  there  is  no  boundary  interference.  Such  may  be  the  case 
in  the  deep  ocean,  where  the  boundaries  are  far  enough  from  the  array  so  that  the 
reflections  off  of  them  are  highly  absorbed.  Also,  if  an  interference  canceling 
scheme  is  used  prior  to  the  Bayes  detector,  the  processor  may  assume  that 
interference  was  canceled  out  completely,  and  implement  the  following  approach: 
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The  processor  observes  a  vector  time  series  as  follows: 


rr(t)  =  [i'0(0.»T(f) . rK_  ,(t)]  - 


-T  T 

-  <  (  <  — 

2  2 


(3.25) 


Under  H, 


rk(t)  =  sk(t)  +  nk{t) 


-T  T 

—  <  t  <  ~ 
2  2 


(3.26) 


Under  H0  : 


rk(t)  =  nk(t)  - 


-T  T 

- <  t  <  — 

2  2 


(3.27) 


As  shown  before,  one  can  write: 


2  1 

nk(t)  =  2  Re  £ak(H)(— )1/2expO2nw0r) 


-T  T 

-  <  t  <  — 

2  2 


(3.28) 


2  1 

sk(t)  =  2  Re  1£bk{n)(—)ir2exp{j2nu0t) 


-T  T 

—  <  r  <  — 
2  2 


(3.29) 


2  i 

rk{()  =  2  Re  ]>>*(«)(— )l/2expO’2nw0t) 


-T  T 

-  <  f  <  _ 

2  2 


(3.30) 


When  noise  alone  is  observed,  Z  .  defined  before,  is  a  zero  mean  complex 
Gaussian  vector  with  covariance  matrix  £. 

When  a  signal  is  present,  Z  is  a  complex  Gaussian  random  vector  with  covariance 
matrix  Q  and  mean: 
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mT  =  [b0{N  Wj) 


bK^(N,).bQ(Nl+  1) . bK_l(Nl  +  1) . ^,,(^2)1  (3.31) 


(3.32) 

(3.33) 


P(Z/HU  •  ,  .  i  , 

A(Z)  = - =  exp[  -(Z-m)  a  (Z  -  m)  +  Z  Q  Z)  (3.34) 

P(Z/H0) 

Any  monruone  function  of  the  Ih'dihood  functional  can  serve  as  a  sufficient 
statistic: 


lnA(Z)  =  -  (Z  -  +  z'a^Z  (3.35) 

=  -  Za~\Z  -m)  +  ma~\z  -  m)  +  z'jCT'z 

=  -  z'a~'z  +  z*a_12H  +  m'or'z  -  m’or'm  +  z\a~!z 

=Z  QTlm  +  m  Q.~lZ  -  m'Q~'m 

since  m  isn’t  related  to  the  observed  vector  Z,  another  sufficient  statistic  is: 

-[Z’jG"‘in  +  m’Q-lZ]  =  Re  [ZO~'m]  (3.36) 

2 

S2 

=  Re  Zz'(n)Q-\n)[b0(n) . bK  _  ,(n)]r 

", 

where  Q~'(n)  is  the  inverse  of  the  nth  block  Q(n)  of  the  block  diagonal  matrix  O- 
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Thus,  the  resultant  optimum  detector  instructs  us  to  implement  a  matched  filter 
where  we  correlate  the  observed  vector  with  the  signal  mean  vector. 

The  above  derivation  serves  as  an  ’example  recipe',  which  is  repeated  in  this 
work  in  order  to  arrive  at  Bayes  optimum  detectors  which  apply  to  various 
scenarios  and  statistical  assumptions. 


4.  JCLSL  and  JCARMA  -  all-zero  and  pole-zero 
adaptive  filters 

One  approach  to  the  problem  of  interference  rejection  in  an  array 
processing  context  is  the  use  of  an  adaptive  beamfor  4ng  structure  which 
dynamically  steers  a  spatial  null  towards  the  source  of  undesirable  signal. 
Interference  leaks  into  the  output  of  a  conventionally  formed  beam  through 
sidelobes  pointing  towards  the  source  of  interference.  One  or  more  reference 
beams  can  be  constructed  which  receive  well  the  interference  but  are  prevented 
from  passing  signals  propagating  from  the  desired  look  direction.  The  outputs  of 
these  reference  beams  are  adaptively  filtered  to  provide  a  good  estimate  of  the 
contaminant,  and  are  then  subtracted  from  the  output  of  the  conventional  beam. 
Such  a  structure  is  performing  a  noise  canceling  operation. 

Since  a  beamformer  weights  and  sums  the  outputs  from  a  finite  number  of 
array  elements,  it  implements  a  spatial  FIR  digital  filter.  The  response 
characteristics  of  a  beamformer  will  exhibit  nulls  as  a  function  of  either:  (1)  source 
arrival  angle  at  a  given  frequency  or  (2)  source  frequency  at  a  given  arrival  angle. 
Beam  based  adaptive  beamformers  can  perform  poorly  due  to  nulls  in  the 
response  characteristics  of  the  reference  channel  beamformer.  A  reference  channel 
adaptive  filter  implementing  poles  can  partially  compensate  for  these  nulls. 

4.1.  An  all-zero  adaptive  filter 

The  all-zero,  joint  process,  least  squares  lattice  structure  (JCLSL)  used  in  this 
dissertation  has  been  discussed  in  detail  elsewhere21’2  along  with  applications  to 
the  adaptive  array  processing  problem.  Figure  4.1  depicts  the  general  structure  of 
this  adaptive  filter,  and  Figure  4.2  details  the  single  stage  structure. 


"3 
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Figure  4.2.  ith  stage  of  the  filter 

The  reference  channel  process  y(n)  is  filtered  to  form  an  estimate  of  the  primary 
channel  process  x{n)  (the  desired  signal).  Of  particular  importance  in  noise 
canceling  applications  is  the  residual  e*M(n)  obtained  by  adding  the  filtered 
reference  channel  to  the  primary  channel.  In  an  adaptive  beamforming  application 
x{n)  is  the  output  of  a  conventional  beam  containing  an  interference 
contamination,  and  y(n)  is  the  output  of  the  formed  reference  beam.  The  filter 
attempts  to  form  a  good  estimate  of  the  contaminant  present  in  x(n). 

4.2.  The  pole-zero  complex  adaptive  joint  process  least  squares  lattice 

This  section  will  formulate  and  solve  the  joint  process  estimation  problem, 
where  the  estimated  process  is  assumed  to  be  an  ARMA  process,  and  therefore  the 
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linear  predictor  takes  the  form  of  a  pole-zero  structure.  The  single  channel  all¬ 
zero  linear  predictor  and  the  all-zero  joint  process  estimator,  as  well  as  the  pole- 
zero  linear  predictor  have  all  been  developed  and  documented2i-2.32, i7.ic.C2_  An 
embedding  approach  is  taken,  where  a  scalar  ARMA  process  is  embedded  into  a 
vector  AR  process,  and  the  previously  developed  AR  recursions  are  utilized  to 
produce  the  pole-zero  joint  process  estimator  recursions.  The  benefit  of  using  the 
more  complex  adaptive  filter  possessing  both  poles  and  zeros  is  finally 
demonstrated  through  a  simulation  example. 

4.3.  The  input-output  ARMA  process  model: 

Consider  a  scalar  data  sequence  y(n)  which  is  assumed  to  be  an  M,h  order 
ARMA  process  of  the  form: 

y(n)+aMAy(n-l)+ . +aMMy(n-M)  =  b„0u(n)+ . +bMMu{n-M)  (4.1) 

where  the  orders  of  the  AR  and  MA  parts  of  the  ARMA  process  are  chosen  to  be 
equal  for  presentation  convenience,  and  u(n)  is  the  input  to  the  ARMA  filter.  Note 
that  by  simply  dividing  equation  (4.1)  by  frM0,  it  can  be  turned  around  to  form: 

u(n)+dM  jU(n-l)+ . +dMMu(n-M)  =  cMOy(n)+ . +cMMy(n-M)  (4.2) 

Here,  u{n)  is  also  represented  as  an  ARMA  process. 

The  innovations  of  yin)  and  u(n)  are,  respectively: 


(”)  =  v(n)-v(n) 

(4.3) 

(UM(n)  =  u(n)-u(n) 

(4.4) 

where 
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y(n)  =  -aMly(n~  1)- . -aMMy(n-M)+bM1u(n- 1)+ . +bMMu(t-M)  (4.5) 

and 

=  ~dM,i  —  1  /  . '  . . M)  (4.6) 

4.4.  The  embedding  approach: 

Equations  (4.3)-(4.6)  can  be  combined  into  a  single  matrix  equation  in  the 
form  of  a  2-dimensional  AR  linear  predictor  equation: 


y(n) 

M 

aM.l 

-K.i 

y(n-i) 

*»• 

= 

u(n) 

/_  i 

■  CM,t 

dM.l  ■ 

.u(n-i). 

Thf  backward  linear  prediction  equation  can  also  be  written  in  terms  of  an 
embedded  ARMA  process  in  a  vector  AR  process: 


y(n— M— 1)1  m 

eM.M+l-U 

~  /m.m+i-u  j[y(n-/) 

-T 

it 

s' 

II 

u(tt_M-i)rs 

i-i 

Vm+1-i  JLu(n-/). 

(4.8) 

where  rY(n)  and  r“(n) 

are  the  Mth 

order  backward  prediction  error 

signals. 

Equations  (4.7)  and  (4.8)  are  then  just  AR  forward  and  backward  linear  prediction 
equations  for  the  2-dimensional  process: 


Y(n) 

u(n) 


(4.9) 


The  most  crucial  assumption  made  here  is  that  u(n)  can  indeed  be  represented  by 
equation  (4.2).  Lee31  shows  that  in  most  cases  this  assumption  is  valid. 
Furthermore,  implicit  in  this  approach  is  the  availability  of  the  innovation  process 
u(n)  to  the  ARMA  filter.  Since  this  input  is  unavailable,  a  bootstrapping  approach 
is  taken,  in  which  <YM,  the  best  estimate  of  u(n),  is  fed  back  to  the  input,  and  used 
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instead.  Now  one  can  use  the  previously  derived  AR  process  recursions  to 
compute  the  ARMA  process  linear  predictor  recursions  and  also  the  joint  ARMA 
process  recursions.  The  ARMA  process  linear  predictor  recursions  are  well 
documented3 1'53  32  and  we  will  not  repeat  them  here.  Instead,  we  will  extend  the 
joint  process  AR  estimator  recursions  to  get  the  joint  ARMA  process  estimator 
recursions.  The  notation  used  in  the  following  equations  is:  y(n)  represents  the 
reference  process  being  filtered  to  predict  the  "desired"  signal  x(n).  u(n)  is  the 
reference  signal  innovation  process.  The  first  index  always  represents  order 
number,  and  the  second  represents  time. 

The  recursion  for  the  desired  signal  prediction  error  in  the  joint  AR  process 
estimator  is31  : 


*  —  *  \*  D~r 

ii+ l.n  — 


(4.10) 


where  is  the  Mth  order  AR  process  cross  correlation  coefficient,  and  R^r  is  the 
Mth  order  AR  process  backward  prediction  error  covariance  inverse.  Using  the 
embedding  approach,  we  now  have  in  the  ARMA  case: 


ry 

ri+l,n 


yi+l,n 


(4.11) 


where,  again,  u(n)  is  the  innovation  process  of  the  reference  signal  v(n).  Now, 
combining  (4.10)  and  (4.1 1)  we  get: 
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It  follows  then  that: 
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Further: 
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Now  define  the  joint  process  reflection  coefficients  as: 
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and  finally  the  desired  signal  prediction  error  recursion  is: 
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The  recursions  for  A*y  and  A*u  are  given  by: 
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yu  yu  £/,n  +  l^i+l.n  +  l 

=  U-o-)Af+1.„+— - 

1  -7f.«+i 


where  (1— a)  is  the  exponential  weighting  factor,  and  -7  is  the  likelihood  variable. 
The  recursions  for  Rryy,  /?ruy,  l?ryu  and  RrulJ  are  given  by: 


And  all  the  above  recursions  yield  the  joint  ARMA  process  estimator  structure 
shown  in  Figure  4.3.  Figure  4.4  details  the  single  stage  structure. 


Figure  4.3.  Pole-zero,  joint  process,  least  squares  lattice  filter 


Figure  4.4.  ith  stage  of  the  structure 


4.5.  Simulations 

Two  sets  of  simulations  will  be  used  to  compare  the  performance  of  the  all¬ 
zero  (JCLSL)  and  the  pole-zero  (JCARMA)  least  square  lattice  adaptive  noise 
canceling  structures. 
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Figure  4.5  describes  how  the  input  data  was  generated  for  the  first 
simulation. 


Main  channel  transfer  function 


Reference  channel  transfer  function 


Figure  4.5.  Data  generation  mode!. 

A  white  noise  time  series  w(n)  has  been  passed  through  two  simple  FIR  digital 
filters.  Hmajn(z),  generating  the  main  channel  data  x(n),  consists  of  a  6-point 
rectangular  window.  Href(z),  generating  the  reference  channel  data  y{n),  consists 
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of  a  4-point  rectangular  window.  Since  the  goal  of  the  adaptive  filter  is  to  invert 
the  characteristics  of  the  reference  channel  transfer  function,  Href(z),  and  to 
model  the  characteristics  of  the  main  channel  transfer  function,  this 

example  presents  a  case  where  the  ability  to  adapt  poles  is  important. 

Figure  4.6  displays  error  power  as  a  function  of  time  at  the  output  of  JCLSL 
and  JCARMA  for  filter  order  p-4.  The  results  show  that  JCARMA  has  a  several  dB 
lower  error  signal  power  than  JCLSL.. 


Figure  4.6.  JCARMA  vs.  JCLSL  error  power  vs.  time  (p-4) 

Additional  insight  is  gained  by  looking  at  the  power  spectra  in  Figure  4.7  (JCLSL), 
and  Figure  4.8  (JCARMA).  The  left  hand  pair  of  panels  are  estimates  of  the  main 
and  reference  channel  input  power  spectra  (white  noise  passed  through  Hmain  and 
Href).  Note  the  dip  in  the  reference  channel  spectrum  which  is  centered  on  the 
first  sidelobe  of  the  main  channel  power  spectrum.  It  is  in  this  particular  region 
that  the  adaptive  filter  can  use  the  ability  to  adapt  poles.  The  right  hand  pair  of 
panels  are  estimates  of  the  adaptive  filter  output  and  error  e(n)  power  spectra.  As 
can  be  seen,  JCARMA  is  able  to  properly  match  the  spectral  characteristics  of  the 
main  channel  where  JCLSL  cannot. 
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Figure  4.9.  JCAKMA  vs.  JCLSL  error  power  vs.  time  (p-8) 


( i«r  IdfO 


of  N^bin  CHomo] 


Figure  4.11.  JCL.CL.  p-8 
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4.6.  REVGEN  simulation 

The  second  simulation  was  set  up  using  a  real  ocean  environment,  and 
REVGEN,  a  REVerberation  GENerator  discussed  in  Chapter  2.  This  is  a 
reverberation  simulation  package  which  implements  directly  the  point  scattering 
model30.  The  returns  from  a  large  number  of  discrete  scatterers,  distributed 
randomly  throughout  the  ocean  volume  and  boundaries,  are  summed  coherently  to 
obtain  the  reverberation  time  series. 

The  sonar  receiving  array  is  made  of  four  elements,  and  is  placed  at  a  depth 
of  150m,  in  water  which  is  450m  deep.  The  array  platform  is  moving  at  12.5  knots. 
A  pure  tone,  180  msec  pulse  is  transmitted  and  propagated  through  the  water 
column.  The  transmit  beampattern  is  uniform  over  60  azimuthal  degrees  and  150 
vertical  degrees.  The  return  volume,  surface  and  bottom  reverberation  time  series 
for  each  of  the  array  elements  is  generated  by  REVGEN.  Two  conventional  beams 
are  formed,  a  primary  beam  using  the  sum  of  all  four  array  element  time  series, 
and  a  reference  beam  using  the  difference  between  the  two  middle  element  time 
series. 

it  is  clear  from  the  above  description,  and  from  examining  the  beampatterns 
that  the  sonar  array  will  receive  boundary  reverberation  which  might  hinder  its 
ability  to  receive  well  from  the  desired  look  direction  (straight  ahead).  An 
adaptive  filter  continuously  steering  a  spatial  null  in  the  direction  of  the 
interference  may  improve  the  sonar’s  performance.  Here,  the  performance  of 
JCARMA  will  be  compared  again  to  that  of  JCLSL,  under  the  described  scenario. 

Figure  4.12  details  the  polar  beampatterns  of  the  main  and  reference 
channels  along  a  vertical  cross-section.  Figure  4.13  details  the  spatial 
beampatterns  of  the  main  and  reference  channels.  Figure  4.14  depicts  range- 
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Doppler  maps  of  the  main  and  reference  channel  time  series,  after  the  platform’s 
own  speed  has  been  compensated  for  by  shifting  the  entire  range-Doppler  map  by 
the  proper  amount.  The  surface  reverberation  return  is  clearly  seen,  peaking 
around  0.5  sec  and  -0.2  cycles/sample. 


Figure  4.14.  Main  and  Reference  channels  range-Doppler  maps 
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Now  we  pass  the  2-beam  vector  time  series  through  the  two  adaptive 
filtering  structures  (JCLSL  and  JCARMA),  and  examine  the  results.  For  both  filters, 
filter  order  was  set  to  3.  and  a  to  0.02.  Figure  4.15  depicts  the  two  adaptive  filter 
outputs.  Though  it  is  clear  that  both  filters  have  done  a  good  job  of  reducing  the 
surface  interference  level,  thereby  allowing  better  reception  from  the  look 
direction,  it  is  not  clear  in  this  case  that  JCARMA  performs  better  than  JCLSL. 
Plotting  the  overall  time  evolving  power  levels  (Figure  4.161  and  the 
JCLSL/JCARMA  error  power  ratio  (Figure  4.17)  further  reveals  that  in  this  scenario 
there  is  only  a  slight  performance  advantage  for  JCARMA  over  JCLSL.  The  error 
power  ratio  departs  0  towards  the  positive  side  wherever  JCARMA  does  better 
than  JCLSL.  The  scenario  used  above  will  be  used  again  in  this  dissertation  to  lest 
and  compare  the  performance  of  different  detectors.  The  above  serves  as  an 
indicator,  which  has  been  confirmed  in  other  simulations,  that  in  these  scenarios, 
the  slight  performance  advantage  of  JCARMA  over  JCLSL,  does  not  justify  the 
higher  algorithmic  complexity,  and  the  higher  run-time  investment. 


PttAJT  4£>) 


5.  Three  specific  problems  and  their  corresponding 
Bayes  optimum  detectors 


This  chapter  deals  with  a  specific  scenario  which  an  active  sonar  system 
may  encounter.  Bayes  optimum  detectors  are  developed  for  three  different  sub¬ 
scenarios,  whose  level  of  uncertainty  is  ever  increasing.  The  setting  is  as  follows: 
an  active  sonar  system  is  mounted  on  an  underwater  platform  which  is  submerged 
in  shallow  ocean  water.  In  this  context  shallow  water  means  that  the  first 
reflection  of  acoustic  energy  back  into  the  receiver  is  significant.  The  ocean  depth 
is  chosen  such  that  there  is  a  significant  no-overlap  zone  between  the  two 
boundary  (surface  and  bottom)  interferences,  and  therefore  only  one  of  them  is 
considered,  meaning  detection  is  performed  in  the  no-overlap  zone.  The 
platform's  depth  in  the  water  column  is  either  known,  or  has  a  known  probability 
density  function.  Sound  speed  profiles  are  assumed  isovelocity,  i.e.  acoustic 
energy  propagates  through  the  medium  in  straight  lines.  The  platform  is 
completely  stationary  in  the  water  column,  i.e.  the  effective  own  Doppler  speed  is 
zero.  The  sonar's  front  end  is  a  four  element  sensor  array  whose  preformed  beam 
outputs  are  available  to  the  processor  in  the  form  of  a  sampled  vector  time  series. 
The  preformed  beam  ensemble  contains  a  sum  beam  which  is  formed  by  summing 
all  element  outputs,  and  one  or  more  difference  beams  formed  by  subtracting  two 
adjacent  element  outputs.  It  is  further  assumed  that  the  transmission  pulse  is 
narrowband,  that  the  incoming  signal  characteristics  are  completely  known,  and 
that  the  target  is  stationary.  The  sonar  is  assumed  to  be  limited  in  performance 
by  acoustic  reverberation  and  not  by  ambient  ocean  noise.  This  assumption  is 
typically  valid  in  medium  and  high  frequency  systems,  certainly  when  dealing  with 
relatively  short  ranges. 

The  three  sub-scenarios  are  as  follows: 
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1.  Platform  is  at  a  fixed,  precisely  known  depth. 

2.  Platform  is  at  an  unknown  depth,  but  the  depth  probability  density  function  is 
known  precisely. 

3.  Platform  is  at  an  unknown  depth,  but  the  depth  probability  density  function  is 
known  precisely.  Also,  there  is  a  fixed  interference  source  whose  direction 
probability  density  function  is  known  precisely. 

Sin^e  the  transmitted  energy  is  propagating  through  the  water,  the 
boundary  interference  direction  is  changing  constantly.  In  the  case  where  platform 
depth  is  known  exactly,  the  interference  direction  is  also  known  exactly.  When 
there  is  uncertainty  in  platform  depth,  it  translates  to  uncertainty  in  the  direction 
of  the  incoming  boundary  interference. 

In  the  first  two  sub-scenarios,  a  sum-difference  beam  arrangement  is 
sufficient,  since  there  are  two  degrees  of  freedom,  namely  signal  and  boundary 
interference.  In  the  third  sub-scenario,  another  interference  source  is  introduced, 
and  therefore  another  degree  of  freedom  is  required.  That  is  why  in  this  scenario 
another  difference  beam  is  formed  and  used,  and  the  problem's  dimension 
increases. 

As  discussed  in  Chapter  1,  ocean  acoustic  reverberation  has  been  the 
subject  of  numerous  studies36-35- !3-46.  "^hese  studies  clearlv  show  that  the 
complex  envelope  magnitude  of  acoustic  reverberation  is  Rayleigh  distributed, 
and  that  the  real  and  imaginary  components  of  this  envelope  are  jointly  Gaussian 
distributed.  This  property  of  reverberation  justifies  using  Gaussian  probability 
density  functions  in  all  the  derivations  of  the  Bayes  optimum  detectors  in  this 
dissertation. 

In  Chapter  3,  a  general  expression  was  derived  for  the  likelihood  function 
in  the  case  of  a  known  signal  in  spatially  correlated  Gaussian  noise: 
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p(Z/H  j) 

A(Z)  = - =  exp[  -  (Z  -  m)'a"‘(Z  -  m)  +  Z  a_1Z]  (5.1) 

p(Z/H0) 

Equation  (5.1)  will  serve  as  the  basis  for  computing  the  Bayes  optimum  detectors 
in  the  three  sub-scenarios  addressed  in  this  chapter.  It  is  clear  that  key  to  this 
derivation  is  our  ability  to  precisely  and  accurately  arrive  at  the  covariance  matrix 
inverse  expressions. 

Figure  5.1  depicts  the  power  spectrum  of  a  typical  sonar  ping  at  the  range 
where  detection  is  attempted. 
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Figure  5.1.  Sum  beam  power  spectrum 

It  is  apparent  that  due  to  this  problem’s  setting  (i.e.  a  narrowband  transmission 
pulse.  0  vehicle  speed,  0  target  speed),  almost  all  the  spectral  content  is  located  in 
frequency  index  0.  or  DC.  For  this  reason,  it  is  sufficient  to  use  a  single  frequency 
term  in  computing  the  received  signal  Fourier  transform  vector  Z.  the  covariance 
matrix  d  which  ends  up  containing  a  single  block,  and  the  signal  mean  vector  m- 
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5.1.  Derivation  of  the  two  beam  covariance  matrix 

Typically,  active  sonar  systems  operating  at  medium  and  high  frequencies, 
are  reverberation  limited.  The  detection  ranges  of  these  systems  are  usually  such 
that  reverberation  dominates  the  ambient  ocean  noise  which  is  therefore 
ignored60.  The  ability  to  precisely  evaluate  the  total  reverberation  (volume, 
surface  and  bottom)  covariance  matrix,  is  the  foundation  to  the  correct  derivation 
of  the  likelihood  functional.  That  will  lead  to  the  design  of  a  trr  'ptima! 
detector  in  the  ''Bayesian"  sense. 

Consider  the  four  element  array  which  we  use.  We  have  available  two  beams: 
sum  and  difference.  The  sum  beam  is  a  sum  of  all  four  channels,  and  the 
difference  beam  is  a  difference  between  the  two  center  channels.  Figure  5.2 

depicts  this  arrangement: 


Figure  5.2.  Sum  and  difference  beams 
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Since  in  this  v  rk.  the  mapping  used  to  reduce  the  received  signal  to  a  finite 
dimensional  space  is  the  Fourier  transform,  the  matrix  sought  here  represents  the 
covariance  between  Fourier  coefficients  of  all  frequency  indices,  and  all  array 
indices.  Flere,  only  a  single  frequency  Fourier  coefficient  is  used,  and  therefore  the 
covariance  matrix  is  block  diagonal,  and  only  one  block  of  it  needs  to  be 
calculated.  Since  the  sonar  is  reverberation  limited,  we  need  consider  only  the 
reverberation  components'  contribution  to  the  covariance  matrix. 

Volume  reverberation  results  from  the  scattering  of  acoustic  energy  by 
marine  life  and  particulate  matter  distributed  in  the  water.  Its  element  to  element 
cross-power  spectrum  depends  heavily  on  the  cross-element  spatial  response 
overlap,  but  by  no  means  is  white.  Array  element  to  element  covariance  is 
therefore  far  from  being  negligible.  It  will  be  shown  here,  though,  that  in  the 
particular  case  of  surn-difference  beams,  volume  reverberation  can  be  assumed  as 
uncorrelated  element  to  another  element. 

Denote  the  k- th  element  received  waveform  by  rk(t). 

Then  the  sum  beam  waveform  5(0  is: 

5(0  =  r0(t)+rx(t)+r2(t)+r3(t)  (5.2) 

and  the  difference  beam  waveform  is: 

d{t)  =  r  ,(t)-r2(t)  (5.3) 

The  cross-covariance  between  the  sum  and  difference  beams  is: 

0^(0  =  E[s(t)d'(t) ]  =  E  [r0(0+M0+'T(0+''3(0](r1(0-r2(0)  (5.4) 

where  £[  ]  represents  the  expected  value  operator,  and  represents  complex 
conjugation. 
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Dropping  the  time  dependence  in  the  notation: 

Qtd  =  E  ror\+riri+r zr\+rir\~>'Q'r'2-r./2-r2r':L-rir\ 


(5.5) 


On  the  average,  and  certainly  when  the  reverberating  volume  is  in  the  far  field  (i.e 
its  range  is  much  larger  than  the  array  length),  it  is  true  that: 


h\ 

=  E 

r2r\ 

«  E 

r*r\ 
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*  E 

rA 

»  E 

r/2 

(5.6) 


and  therefore: 


Q*d 


(5.7) 


We  see  that  in  the  particular  case  of  sum-difference  beams,  far  field  volume 
reverberation  is  essentially  uncorrelated,  and  its  only  contribution '  to  the 
covariance  matrix  is  along  the  main  diagonal. 

Consider  now  the  boundary  reverberation.  Denote  the  nth  Fourier  coefficient 
at  the  0th  array  element  due  to  the  boundary  reverberation  as  d0{n).  Then  the 
Fourier  coefficient  at  the  pth  array  element  is: 


dp(n)  =  d0(n)exp\-jnu06p}  P= 0....3  (5.8) 

where  Sp  is  the  time  delay  between  element  0  and  element  p.  The  sum  Fourier 
coefficient  is  then: 

3 

ds(n)  =  £  d0(n)exp\-jn~-05p}  (5-9) 

P  -  0 


and  the  difference  Fourier  coefficient  is: 


dd(n)  =  d0(»)(exp(-J/f-.-0$,)-exp(-./M; 


0^1 


Consider  now  the  covariance  matrix,  let  exp(jnu06x)  =  W.  and  note  that 
then  the  cross  covariance  is: 


Q„,  =  E\ds(n)dd(n)}  =  |  dQ(n)  |  ^l  +  V^'+W^+vr'jfw'-W2) 


=  I  d0(n ) 


1  jvr2-^2) 


=  -2j  |  d0l  ‘S!n(2«w05,) 


In  general,  if  there  are  k  (even)  elements  in  the  array,  and  the  middle 
subtracted  to  form  the  difference  beam: 


E\ds{n)d'd{n)]=  -2j\  d0  \  2sin(m^0— 6,) 
t  J  2 

x-i  k- 1 

When  k  is  odd,  and  elements  subtracted  are - and - +1: 
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The  sum  beam  autocovariance  is: 


Qss  =  £[dJ(n)d>)]=  Y.doWVT’Y.doWVf 

p— o  p^> 


=  |  d0 1  2(w3 +2W  ->-3Wl +4+3  W~l+2W~2+W~3 


(5.10) 
&p  =  PS  i. 

(5.11) 


two  are 

(5.12) 

(5.13) 

(5.14) 
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=  2  |  d0\  2|4+3cos(/iw0<51)4-2cos(2m'Jj0<5,)+cos(3«^0<51)| 

And  in  general: 

jfc-i 

E  ds(n)d'(n) j=  2j  d0  \  2 '£(k-i)cos(inu06i) 

1-0 

The  difference  beam  autocovariance  is: 

Qm  =  £[^(«)^'(«)]=  i  d0\  2  (( W1- W^2)( W1- W2))  (5.15) 

=  2  |  d0  \  2[l-cos(Mw051)j 

The  covariance  matrix  can  now  be  inverted  and  used  in  the  likelihood  ratio 
expression. 

S.2.  Signal  known  exactly  in  correlated  Gaussian  interference  of  precisely 
known  direction 

The  covariance  matrix  elements  computed  are  all  functions  of  the  basic 
element  to  element  time  delay.  When  depth  is  known  precisely,  figuring  out  <5, 
requires  a  trigonometric  computation  which  is  illustrated  in  Figure  5.3: 


Ocean  surface 


Figure  5.3.  Element  to  element  time  delay  computation 


sina  = -  (5.16) 

Z(T) 

Z(T)  =  cT 
d 

sina  =  — 
cT 


and  therefore: 


x(5t)  =  (sina 


*(*i)  /  sina  Id 


(5.17) 


(5.18) 


where  /  is  the  inter-element  spacing,  d  is  the  array  depth,  c  is  the  speed  of  sound, 
and  T  is  the  observation  time. 


The  following  summarizes  our  findings  relative  to  a  representative  block  of 


the  block  diagonal  reverberation  covariance  matrix: 
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k-1 

2  I  dQ\  2'£l(k-i)cos(n<j0i6l)+Vs 

i-0 


i  i  2  * 

-2j  j  d0  I  sin(nu0~6,) 
7 


Qin.Sj)  = 


2  j  I  d0 1  2sin(« 


Vj+2  |  d0|  2  |l-cos(«'^06,) 


(5.19) 


Whcic  Vs  and  Vu  are  the  sum  and  difference  volume  reverberation  powers. 

Now  one  can  invert  the  covariance  matrix,  and  use  it  in  the  previously 
computed  likelihood  ratio  expression  derived  for  this  case  of  signal  known 
exactly  in  correlated  Gaussian  interference  of  known  location. 


piling 

A (Z)  = - =  exp[  -  (Z  -  m)  Ci\Z  -  m)  +  Z'a_12]  (5.20) 

P(Z/H0) 

5.3.  Signal  known  exactly  in  correlated  Gaussian  interference  of  uncertain 
direction 

Boundary  interference  in  this  case  is  coming  from  an  uncertain  direction 
since  vehicle  depth  is  uncertain.  We  do  have  some  limited  knowledge  about  depth, 
though,  in  the  form  of  a  probability  density  function  (pdf)  of  it.  This  pdf  is 
translated  directly  (through  a  simple  transformation)  to  a  pdf  of  the  element  to 
element  delay  <5j.  The  same  covariance  matrix  expressions  developed  for  the 
previous  case  still  apply,  but  now  computing  the  likelihood  ratio  involves 
integration  over  the  uncertain  parameter.  When  H0  is  in  force,  the  conditional 
probability  density  function  is: 

p{Z/H0,&x)  =  (xk  |  Q(n.5,)|  r'expj— Z'(n)a_1(n,i,)Z(n))  (5.21) 

In  our  case  Z  is  2-dimensional  (sum-difference),  and  only  one  frequency  index  is 


used.  Therefore: 
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Z  = 


(5.22) 


The  covariance  matrix,  which  is  a  function  of  i5,  has  been  derived  previously.  Let: 


then: 


Qiti) 


Q.II  Ql2 

Qzi  Ql2 


(5.23) 


Eo  =  —  Z  Q  (5[)Z  =  —  ( I  2 1 1  Qz2+  I  2 2 1  Qj  j— 2Re(z2z i Q^i)]  (5.24) 

A  v  ' 

where  A  is  the  determinant  of  Q.  Note  that  A  ,  and  Q,,  are  all  functions  of  <5,.  Now 
the  conditional  pdf  under  H0  assuming  5,  is  known  can  be  written  as: 


p(Z/H0.&x)  =  ~~exp(E0)  (5.25) 

r  A 

The  knowledge  about  the  distribution  of  the  element  to  element  delay  can  now  be 
utilized  to  form  the  pdf  under  H0: 


p(Z/H0)  =  — -  f  exp(E0)p(Sl)Z51  (5.26) 

*  A  «„ 

where  Su  and  8lf  represent  the  boundaries  in  between  which  p(<5,)  is  defined. 
Finally,  when  assuming  uniform  distribution  on  5,,  as  will  be  assumed  here: 


P(Z/H0)  = 


— - /  exp  (E0)d8x 

x  A(<5(  f—8u)  ^ 


(5.27) 


When  Hj  is  in  force,  the  conditional  probability  density  function  is: 
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P(Z/H =  (/  |  Q(n,5x)  |  )  'exp  (-  fz(n)-60(n)a(n)] 

^2_'(«.<5i)  [Z(n)-b0(n)u(n)^ 

where  b0(n)  is  the  signal  mean  vector,  and  u(n)  is  the  signal  "pointing 
signal  direction  is  assumed  precisely  forwards,  i.e.: 


u  = 


Using  the  same  steps  used  to  derive  p{Z/HQ,8x)  : 

[-(|  zt !  2Q22+i  Z2I  2Qn-2Re(z;zlQ21)) 


_1_ 

A 


+2fa0Re(z;Q22-z;Q21-z;Q12+z;Qu)-|  b0 1 


The  conditional  pdf  under  H x  assuming  5,  is  known  can  be  written  now 


p{Z/H ,  <5j)  =  -7-exp  (E,) 
/a 


Using  what's  known  about  the  distribution  of  6,: 


°\r 


p(Z/Hx)  =  — —  J  exp(El)p(6l)d5l 


?r  A 


and  when  the  uniform  distribution  is  incorporated: 


° if 


p(Z/H ,)  = 


jt  A(<5[  f-S ,,) 


Jexp(E1)dSl 


(5.28) 


vector.  The 


(5.29) 


(5.30) 


>1 

as: 

(5.31) 


(5.32) 


(5.33) 


We  are  now  finally  ready  to  compute  the  likelihood  ratio  for  the  case  of  signal 
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known  exactly  in  correlated  Gaussian  interference  of  uncertain  direction: 


/  exp  {El)d6l 
P(Z/HQ  ~  ^ _ 

P(Z/H0 )  s\f 

/  exp(E0)dSl 


(5.34) 


S.4.  Derivation  of  the  three  beam  covariance  matrix 

When  we  consider  a  more  complex  case,  in  which,  in  addition  to  the 
boundary  reverberation  interference  we  have  another  fixed  interference  source 
(e.g.  a  jammer),  we  need  another  degree  of  freedom  in  the  system.  Therefore, 
another  difference  beam  is  formed,  and  we  have  three  beams:  one  sum  beam  and 
two  difference  beams.  Figure  5.4  depicts  this  arrangement: 


Figure  5.4.  Sum  and  difference  beams 
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Now  we  can  attempt  to  compute  the  increased  order  covariance  matrix.  Denote 

r/f  th 

the  n  Fourier  coefficient  at  the  0  array  element  due  to  the  boundary 

th 

reverberation  as  d0(n).  Then  the  Fourier  coefficient  at  the  p  array  element  is: 

dp(n)  =  d0(n)exp|-jnu;0«5pj  P=0,..,3  (5.35) 

where  $p  is  the  time  delay  between  element  0  and  element  p.  The  sum  Fourier 
coefficient  is  then: 

3 

ds(n)=  £  d0(n)exp[-jnu06^  (5.36) 

P  -  o 

and  the  difference  Fourier  coefficients  are: 


ddl(n)  =  (d0(«)[exp|-jnw050j-exp[-jntj0<51jj  (5.37) 

dd2{n)  =  d0(n)  (exp  (-j/iw0<52J-exp  (-jnw0<53]}  (5.38) 

Consider  now  the  covariance  matrix,  let  exp(jinu;0i51)  =  W,  and  note  that  6p  =  pSv 
then  the  first  cross-covariance  is: 

=  E[ds(n)ddl(n)]  =  |  d0{n) !  +  (5.39) 

=  |  cf0(n)  |  2(w^3-Wl) 

The  second  cross-covariance  is: 

0jd2  =  6 [ds(n)dd2(n)]  =  |  d0(n) !  2[\+vr\vr2+vrz){wi-wi)  (5.40) 

=  I  dQ(n)\  2[ur'-w3) 
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The  sum  be,  itocovariance  (calculated  previously)  is: 

*-i 

Qts  =  E  ds{n)d](n)  =  2|  d0  \  2  ^(k-ficosUn^J  (5.41) 

The  difference  beam  autocovariances  (calculated  previously)  are: 

Q^ui  =  E  ddl{n)d'dl{n)  =  E^dd2{n)dj2(n)  =  2  |  d0  \  2  [l-cos(n^1)j  (5.42) 
The  difference  beam  cross  covariances  is. 

I  dQ(n)\  2(l-W^,)(lV2-IV3)  (5.43) 

=  |  d0(n)  \  2(2W*-Wl-W3) 

When  a  fixed  interference  is  introduced,  the  structure  of  the  resulting  covariance 
matrix  is  identical  to  the  boundary  reverberation  covariance  matrix.  The  only 
difference  is  that  the  time  delay  variable  6  does  not  vary  with  range,  but  is  fixed. 
The  fixed  time  delay  is  denoted  5f.  The  combined  covariance  matrix  is  then  the 
sum  of  the  range  variable  boundary  interference  matrix,  the  fixed  interference 
matrix,  and  the  diagonal  volume  reverberation  matrix: 

On  0-12  On 

Q.  —  f)  =  Qz i  Qi 2  O23  (5.44) 

■Q31  O32  Q.33 

5.5.  Signal  known  exactly  in  correlated  Gaussian  interference  of  uncertain 
direction,  plus  a  fixed  interference  of  uncertain  direction 

Fixed  interference,  as  well  as  boundary  reverberation  in  this  case  are  coming 
from  uncertain  directions.  Again,  we  do  have  some  limited  knowledge  about  fixed 
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.nterference  direction  in  the  form  of  a  probability  density  function  (pdf)  of  it. 
Here,  computing  the  likelihood  functional  involves  integration  over  both  uncertain 
parameters,  namely  and  8f.  When  H0  is  in  force,  the  conditional  probability 
density  function  is: 

p{Z/H0,b\,bf)  =  (/  |  Q.(n,5lt6f)\  r‘exp[-Z  (n)Q~1(n,8l,8  r)Z(n)j  (5.45) 

In  our  case  Z  is  3-dimensional,  and  only  one  frequency  index  is  used.  Therefore: 

Z  =  z2  (5.46) 

,z3J 

The  covariance  matrix,  which  is  a  function  of  8t  has  been  derived  previously.  Let: 

Hit  Hj2  H13 

Q~'(S,,Sr)  =  =  H2I  H22  H23  (5.47) 

h3.  h32  h33. 

then: 

f0  =  -ZH(8vSf)Z  (5.48) 

=  -(|  z,  |  2Hn+i  z2|  2H22+|  Zj|  2H33+2Re(z2zlH2l)+2Re(z'3zlH3l)-i-2Re(z'3z2H32)} 

Note  that  H(j  are  functions  of  <5,  and  Sf.  Now  the  conditional  pdf  under  H0  when 
and  8 f  are  known  can  be  written  as: 

1 

p(Z/H0,i]8  f)  =  exp(£0)  (5.49) 


Where  A  is  the  determinant  of  Q,  and: 
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A  =  MSlmSf)  (5.50) 

Co  =  Co  f) 

The  knowledge  about  th<"  distribution  of  5,  and  6f  can  now  be  utilized  to  form  the 
pdf  under  /  u.  Assuming  that  the  two  time  delays  are  independent,  which  is 
reasonable: 


p(Z/H0 )  = 


1 


T3A  s  : 


J  J  exp(£0)p(5,)p(5/:)d51rf<5/ 


(5.51) 


where  Su  and  5xf  represent  the  boundaries  in  between  which  p(Sx)  is  defined,  and 
if,  and  iff  represent  the  boundaries  in  between  which  p(5 f)  is  defined.  Finally, 
when  assuming  that  both  delays  are  uniformly  distributed,  as  will  be  assumed 
here  we  get: 


P(Z/H0)  = 


When  H]  is  in  effect,  the  probability  density  function  is: 


/  /  exp(E0)dild6f 


(5.52) 


p{Z/Hx,ix,bf)  =  {/ A)  'exp  (—  \Z(n)-b0{n)u{n)^  (5.53) 

Q~'(n,Sl,Sf)  [z(n)-b0(n)li(n)j] 

where  b0(n)  is  the  signal  mean  vector,  and  u(n)  is  the  signal  "pointing"  vec  or.  The 
signal  direction  is  assumed  precisely  forwards,  i.e.: 


U  = 


T 

1 

1 


Using  the  same  steps  used  to  derive  p(Z/HQ,Sx5 f)  : 


(5.54) 
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E\  =  —  ^Z— fa(>ii|  £2  (^i^f)  [z~  (5.35) 

=  —  Eq+2 b^Re  i(H n+h  12+EI 12)+ Z2(H2\J<-H 22+ H 2j)+z 3(H 3i+H 22+ H ji) j 

3  3 

-I  &o 1 2ES^j 

i-lj-l 

The  conditional  pdf  under  Hl  when  both  delays  are  known  can  be  written  now  as: 


1 

p{Z/Hl8[5f)  =  — —  exp(£,) 
ir  A 

Using  what  is  known  .  oout  the  distribution  of  S{. 


P(Z/H{) 


hi*  ft 


x3A 


/  /  txp(E \)p{6\)p(& f)d6xd8 . 


6u6r, 


and  when  the  uniform  distributions  are  incorporated: 


(5.56) 


(5.57) 


PU/H 1)  =  — - /  /  exp {Ex)d5xd5f  (5.58) 

K  M6lf-Su)(8ff-Sfi) 

We  are  now  finally  ready  to  compute  the  likelihood  ratio  for  the  case  of  signal 
known  exactly  in  correlated  Gaussian  interference  of  uncertain  direction  plus  an 
interference  of  uncertain  direction: 


/  /  exp (Ex)d5xd8f 
P(Z/HQ  »»*n _ 

P(Z/H0)  ~  Vrr 

/  /  exp(E0)c/<5,d<5f 


(5.59) 
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We  have  developed  in  this  chapter  three  different  Bayes  optimum  detectors 
for  three  different  ocean  environments.  These  will  be  used  later  in  this 
dissertation  in  the  detector  comparison*:. 


6.  Block  processing  and  time  sequential  approaches 
to  Bayes  optimum  detection 

The  general  expressions  developed  in  the  previous  chapter  for  the  three 
different  cases  are  indirectly  all  time  dependent.  As  the  transmission  pulse 
propagates  through  the  water,  it  ensonifies  different  regions  of  the  ocean,  and  the 
return  signal  varies  therefore  with  time.  The  time  dependence  enters  the  optimum 
detector  expressions  through  the  element  to  element  delay  <5,,  and  through  the 
signal  and  boundary  reverberation  powers. 

There  are  at  least  two  approaches  to  processing  the  data.  The  simpler 
approach  is  to  divide  the  incoming  data  stream  into  equal  length  blocks,  and  to 
process  each  block  separately,  without  carrying  over  any  information  from  the 
processing  of  one  block  to  the  next.  This  approach  is  called  block  processing,  and 
is  obviously  sub-optimal  in  a  sense,  since  it  takes  a  look  at  only  part  of  the 
available  information  at  any  given  point  in  time.  Despite  its  apparent  flaw,  this 
approach  is  many  times  easier  to  implement,  and  is  well  suited  for  real  time 
applications.  The  second  approach  is  to  store  all  the  incoming  data,  and  then 
process  it  at  once  and  make  a  decision.  This  approach  is  called  one-shot,  and  it  is 
optimal  in  the  sense  that  it  makes  use  of  all  the  available  data.  A  variant  of  this 
approach  is  the  time  sequential  approach.  Here,  the  processor  operates  on  a  block 
of  observations  at  a  time,  but  through  updating  the  a-priori  probability  density 
functions  of  the  uncertain  parameters,  it  carries  over  acquired  information  from 
block  to  block.  This  approach  is  mathematically  equivalent  to  the  one-shot 
approach,  and  produces  exactly  the  same  decisions,  but  in  some  respects  is  more 
attractive  since  it  demonstrates  the  inherent  learning  nature  of  the  detector,  and 
lends  itself  to  comparisons  with  other  adaptive  structures.  Both  the  block 
processing  and  the  time  sequential  approaches  have  been  implemented  in  this 
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dissertation,  and  are  both  used  in  the  comparisons  with  ad-hoc  detectors 

6.1.  The  time  sequential  processor  design  equations 

Consider  the  vector  rT(t)  of  time  waveforms  observed  by  the  receiver  array 
as  in  (3.20).  £T(t)  is  broken  into  time  sequences  of  length  Tinc,  and  the  vectors  of 
Fourier  coefficients  zj{n)  are  formed,  where  i  represents  the  i th  increment.  Thus, 
for  a  total  observation  period  of  T  =  LTinc,  t{t)  will  be  represented  by  L  vectors  of 
Fourier  coefficients,  iZ\ . 2L). 

When  using  the  block  processing  approach,  each  one  of  vectors  zj 
represents  a  block,  and  is  being  processed  separately  from  all  the  others.  When 
using  time  sequential  processing,  the  approach  is  as  follows.  In  order  to  form  the 
likelihood  ratio,  we  need  the  marginal  distributions  of  the  observables, 
conditioned  on  the  hypothesis  in  force.  Here,  the  observables  are  the  L  vectors 
(Z . . 2l)-  Suspending  the  conditioning  on  Hl  and  H0  for  the  moment: 

L 

P(Z  . . 2l)  —  TLP(Zi/Zi-v . •Zi)  (6.1) 

i-i 

L 

=  I  n>(Z,/Z,_, . 

i  t- 1 

where  £  is  the  vector  of  uncertain  parameters,  which  for  simplicity  is  assumed 
time-invariant  here.  Assuming  parameter  conditional  independence  of  the  vectors 

Z,: 


p(Z,/Z,_, . =  PiZj ft 


(6.2) 


Substituting  now  (6.2)  into  (6.1)  we  get: 
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P<Z, A)  =  I II Ptfi/&P<Qd&  (6.3) 

i  I- 1 

Now,  applying  Bayes’  rule  L  times  to  the  integrand  in  (6.3): 

L 

p(Zx . A)  =  FI  /  PA/QPtilZi-i . Zx)dL  (6.4) 

i- 1  i 

where  pi^Z^ . Zt)  is  an  updated  version  of  the  a-priori  probability  density 

function  of  £ : 

p(Zi^\, . i>£) 

PWM . A)  -  (6.5) 

p(Z,-i . A> 

p(Zj_i/Zi_2< . uZi^}p(Zj_2> . %ZxJl} 

P<A-\ . -2]) 

and  assuming  parameter  conditional  independence: 


p(Zi_i/0p(Z,_2 . A<& 


pA-x . A) 


(6.6) 


p(Zi_\/fi)p(Ji/Zj-2 . >Z\)p{Zi_2' . Ax) 

p(Z(_i/Zj_2’ . ‘Zi)P(Zj_2 . mZ[) 

p(Zj_i/^)p{SJZj^2 . >Z  i) 

P(Zi-i/Zt-2>“—’Z \) 

Note  that  the  updated  version  of  the  a-priori  probability  density  function  of  £  is 
formed  using  probability  density  functions  which  are  either  known,  or  computed 
in  the  previous  iteration. 


Equations  (6.1),  (6.4)  and  (6.6),  when  conditioned  on  H,  and  H0,  are  the 
design  equations  used  to  obtain  the  marginal  distributions  in  the  likelihood  ratio 
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expression  for  L  iterations: 


P(Z  i Zx J  H  i) 

P(Z  i< . Zi_/H  o) 


L 

n  I  pUJLHJpWZ^ . ZxHx)dk 

/-I  i _ 

L 

n  I  p(Zi/&.H0)pti/Zi_l . ZxH0)di 

i-i  i 


(6.7) 


Note  that  since  we  are  dealing  with  an  active  sonar  problem,  the  signal  (target 
echo)  is  present  only  in  the  final  iteration  currently  considered.  The  assumption  is 
made  that  there  is  only  one  target  in  the  look  direction  and  that  it  resides  in  the 
current  range  cell.  Thus  the  updated  a-priori  probability  density  functions  under 
both  hypotheses  are  identical  in  all  iterations  but  the  last. 

6.2.  Time  varying  uncertain  parameter  vector 

The  previous  discussions  limited  the  uncertain  parameter  vector  to  be  time- 
invariant.  In  actuality,  the  uncertain  parameter  may  be  time-varying.  In  the  first 
case  we  are  dealing  with,  there  is  a  single  uncertain  parameter  -  boundary 
reverberation  element  to  element  time  delay,  which  certainly  is  time-varying.  That 
results  in  having  to  find  a  transition  equation  which  makes  use  of  the  problem 

geometry,  to  update  p{8JZ . ZxHj)  for  iteration  /+1.  The  likelihood  ratio 

expression  is  now: 


PI  J  p{Z(/8j  H i)p(t>i/Z(_{. . Z\  H i)d5i 

pUi . i-l 

=  =  " 

p(Zx . ZJHJ  L  , 

n  I  PiZ'/S'HJp^/Z^ . Zl,H0)dSi 

i-l  s, 

and  the  a-priori  probability  density  function  update  equation  is: 


(6.8) 


p(<5j/Z;_i . Z\,H j) 


p(Zf_i/5/_i.H  j)p(&i/  Z,_2 . Z \'H j) 

pUi.x/Z^ . Z\,Hj) 


(6.9) 
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where  j  =  0,1  represents  the  hypothesis  in  force.  The  need  for  a  transition 
equation  arises  since  following  iteration  i-1.  we  have  available  only 
p(bi-i/Zi_2 . and  in  (6.9)  we  need  p(Sj/Zi_2 . . 

The  transition  equation  is  developed  in  Appendix  B,  based  on  the  geometry  of 
Figure  5.3,  and  the  result  is  quoted  here: 


=  — : -  (6.io) 

l—csl 

1+ - 

/(/- 1) 

Making  use  of  this  transition  equation,  we  can  arrive  at  the  desired  pdf  by 
computing  the  pdf  of  a  function  of  a  random  variable: 


p(<5,_i/Z(_2> Zll,H ,) 

P(W_2 . Zx,Hj)  = - - -  (6.1 1) 

aS, 

'  ds,-t  ' 5"' 

Appendix  B  also  develops  the  exact  expression  for  the  above  derivative. 
Equations  (6.8),  (6.9)  and  (6.11)  constitute  the  final  design  equations  of  the  time 
sequential  Bayes  optimum  processor  with  a  single  uncertain  parameter. 


We  now  examine  the  more  complex  case  where  a  fixed  interference  source 
of  uncertain  direction  is  added  to  the  scenario.  Here,  the  uncertain  parameters 
indeed  form  a  two-dimensional  vector: 


(6.12) 


where  <5,  is  the  time-varying  boundary  reverberation  element  to  element  time 
delay,  and  sf  is  the  time-invariant  fixed  interference  element  to  element  time 
delay.  We  again  have: 
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n/  P{Zi/&.Hl)p(i/Zl_l . ZuHx)dL 

P(Z  1 . ZL/H  i)  /-I  i 

At- - =  - -  (6.13) 

P(Zl ZJH0)  ^ 

n  /  P(Z,/^H0)p(A/Z,_l . ZxH0)di 

/-l  i 

and  since  6t  and  5F  are  independent  random  variables: 

L 

n/ 1  P(Zj/&i  6  f  H  i)p(£j/Zi_i . ’tZi'H  . uZi  H  i)d6jd6  f 

i-1  ss 

\L  =  - -  (6.14) 

FI  //  P(Zl/6i  df  H0)p(Sl/Zl,l, . .Zj ,H0)p(8j/Zj_i . vZi  H0)dSjd5  f 

iml  st6f 


Obtaining  the  update  equations  of  6,  and  8f  is  a  simple  extension  of  the  ideas 
previously  described  and  will  not  be  repeated  here.  Suffice  it  to  say  that  the 
transition  equation  for  <5,  remains  the  same,  while  no  transition  equation  is 
required  for  5f  since  it  is  time-invariant. 
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7.  Detector  evaluation  and  comparison  -  methods 
and  tools. 


7.1.  Detector  comparison 

As  described  in  Chapter  3,  the  performance  of  two  detection  schemes  will 
be  evaluated  and  compared.  One  scheme,  which  represents  an  ad-hoc  approach, 
employs  an  adaptive  filter  to  first  reduce  boundary  interference  level.  It  then 
assumes  that  the  processed  signal  contains  only  the  volume  reverberation 
component,  and  treats  the  problem  as  one  of  signal  known  exactly  in  noise  known 
exactly.  The  adaptive  filter  chosen  to  be  used  in  the  comparison  is  the  all-zero 
least  squares  lattice  filter  (JCLSL).  It  has  been  shown  previously  (Chapter  4),  that  it 
is  quite  adequate  for  the  task,  and  that  a  pole-zero  adaptive  filter  has  only  a  very 
slight  advantage  over  it,  in  the  scenarios  used.  The  second  scheme  is  the  result  of 
a  more  global  approach  which  uses  our  knowledge  of  the  problem's  statistics  to 
arrive  at  an  optimum  solution  in  the  Bayesian  sense.  Figure  7.1  depicts  the  two 
schemes. 
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Figure  7.1  The  two  detection  schemes 

The  array  used  in  the  comparison  is  a  4-element  array,  where  the  main 
beam  was  derived  by  summing  all  four  element  outputs,  and  the  reference  beam(s) 
was  generated  by  subtracting  two  adjacent  element  outputs. 

7.2.  Adaptation  parameter  selection 

One  of  the  parameters  defining  the  behavior  of  an  adaptive  structure  is 
the  adaptation  coefficient  a,  which  takes  values  between  0  and  1.  This  parameter 
determines  the  filter’s  "forgetting  factor",  or  how  much  emphasis  the  filter  puts  on 
past  samples.  Typically,  the  filter’s  time  constant  (in  samples)  will  be: 

1 


Oe 


(7.1) 


When  using  an  adaptive  structure  such  as  the  adaptive  filter  used  in  this 
dissertation,  the  questions  always  arise  as  to  how  is  the  adaptation  coefficient  a 
selected,  and  for  how  long  the  adaptive  filter  should  be  permitted  to  run  and 
adapt  before  its  output  is  considered  valid  to  sample.  In  order  to  resolve  these 
questions,  a  run  was  generated,  where  a  simulated  target  echo  in  the  look 
direction  was  superimposed  on  an  ordinary  reverberation  return  ping.  This  ping 
was  picked  out  of  the  ensemble  of  pings  which  are  later  used  in  this  work  to 
compare  the  different  detection  approaches.  The  ping  was  processed  three  times 
by  the  same  adaptive  filter  with  three  different  adaptation  coefficients,  namely 
0.02,  0.05  ,  and  0.2.  Then,  the  adaptive  filter  output  power  was  plotted  against 
main  beam  power,  for  the  three  different  a’s.  Figure  7.2  details  the  results. 
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MoJn  boom  \im  ■  JC1_SL  outpu-fc  (Plpha  ~  0.02) 


Main  boom  ua,  JC3_St_  oulpul  (Plpha  *  0.05) 


Main  beam  v*  JC1_SL_  output  (Plpha  =  0.2) 


Figure  7.2  Main  beam  vs.  adaptive  filter  output  power. 

It  is  clear  from  observing  the  plots,  that  as  scon  as  a  is  chosen  larger  than  0.02, 
the  filter  starts  canceling  the  target  echo  as  well  as  the  reverberation  interference. 
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This  is  undesirable.  Therefore,  throughout  this  work  a  was  kept  at  0.02.  The 
selection  of  a  answers  the  second  question  as  well,  since  an  adaptive  filter 
typically  takes  a  few  time  constants  to  adapt.  With  this  selection  of  a,  the  time 
constant  r  is  50  samples.  Therefore,  sampling  the  filter’s  output  anywhere  beyond, 
say,  500  samples  (10  time  constants)  is  safe. 

7.3.  Test  scenario  description 

The  two  detection  approaches  are  tested  and  compared  under  three 
different  interference  conditions: 

Case  1.  Boundary  interference  is  coming  from  a  known  direction. 

Case  2.  Boundary  interference  is  coming  from  an  uknown  direction,  whose 
probability  density  function  is  known. 

Case  3.  Boundary  interference  is  coming  from  an  uknown  direction,  whose 
probability  density  function  is  known,  and  there  is  another  fixed  source  of 
interference  coming  from  an  unknown  direction  whose  probability  density 
function  is  known. 

One  source  of  uncertainty  in  boundary  interference  direction  is 
uncertainty  in  the  array  platform’s  depth.  Therefore,  depth  was  chosen  as  the 
uncertain  parameter,  and  its  probability  density  functions  were  assumed  known. 
When  dealing  with  the  fixed  interferer,  probability  densities  of  direction  were 
defined  directly.  In  all  cases  of  uncertainty,  uniform  probability  density  functions 
were  assumed.  The  basic  test  scenario  selected  for  Case  1  above  is  described  in 


Figure  7.3. 
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Figure  7.3  Basic  test  scenario  (Case  1) 

As  can  be  seen  from  the  figure,  the  bottom  depth  is  purposefully  given  a  large 
value,  so  that  bottom  reverberation  will  come  in  late,  and  only  one  source  of 
boundary  reverberation  (surface)  is  dealt  with. 

Figures  7.4  and  7.5  describe  the  test  scenarios  for  Cases  2  and  3  above. 
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Figure  7.5  Test  scenario  for  uncertain  depth  and  fixed  interference 
(Case  3) 

The  output  of  both  detectors  s  examined  approximately  0.2  sec  following 
transmit,  or  at  about  the  800^  sample.  This  gives  the  adaptive  filter  enough  time 
to  adapt.  This  time  corresponds  to  a  distance  to  the  illuminated  surface  patch  of 
150m.  The  maximum  depth  uncertainty  used  in  this  work  is  54m,  centered  around 
30m.  That  spans  a  depth  between  3  and  57m.  The  corresponding  angular  span  to 
the  surface  is  as  follows: 


-i 


<Xrr,.,  =  sin 


54 


I  150 


_-.L! 

150  J 


=  21.1 


=  1.15  s 


(7.2) 


(7.3) 


Similarly,  the  fixed  interference  direction  uncertainty  spans  the  range  of -10°  to 
-20°.  Examining  the  sum  and  difference  polar  beampatterns,  we  see  that  within 
the  above  angular  ranges,  both  are  well  behaved,  and  the  main  beam  has  no  nulls. 
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Figure  7.6  Sum  and  Difference  polar  beampatterns 

7.4.  ROC  curves 

The  output  produced  by  both  detection  schemes  is  the  likelihood  ratio 
which  in  turn  is  compared  to  a  threshold  in  order  to  make  the  decision  whether  a 
target  is  present  or  not.  The  likelihood  ratio  A,  summarizes  all  the  information 
about  the  data  for  a  given  observation  interval  and  therefore  is  significant  in 
evaluating  the  detector's  performance.  A  complete  description  of  a  detector 
includes  the  likelihood  ratio  A  and  the  receiver  operating  characteristic  (ROC) 
curves.  These  curves  plot  the  detection  probability  versus  the  false  alarm 
probability,  where  the  threshold  q  serves  as  a  parameter.  The  false  alarm 
probability  PF  and  the  detection  probability  PD  are: 

OO 

PF  =  Jp(A/H0)dA 


(7.4) 


(7.5) 


PD  =  Jp(A/H,)dA 
n 

Where  r?  is  the  detection  threshold,  it  has  been  shown48  that: 

P(A/H,)  =  \p{\/H0) 

Thus  PD  can  be  written  as: 


(7.6) 


PD  =  jAp(K/H0)dA  (7.7) 

n 

When  the  densities  of  A  cannot  be  determined  analytically,  one  carries  out  a  Monte 
Carlo  simulation  of  the  detector  under  test.  From  the  simulation  results  one  then 
forms  estimates  p{A/Hx)  and  p(A/H0)  of  the  desired  densities.  These  can  be  used 
to  calculate  PF  and  PD  which  constitute  the  ROC  curve.  We  see  that  (7.47  above 
implies  that  only  the  density  under  H0  need  be  obtained.  This  eliminates  the  need 
to  simulate  the  signal.  Note  that  (7.7)  places  an  emphasis  on  the  upper  tail  of 
p(A/H0),  where  typically  only  a  few  observations  will  lie.  An  equivalent  expression 
to  (7.7)  is: 


PD  =  1  -  /Ap(A/H0)dA  (7.8) 

o 

and  it  is  shown  in27  that  it  is  actually  preferable  to  use  (7.8). 

As  a  point  of  reference,  and  an  example  of  ROC  curve  generation,  100 
observation  intervals  of  white  Gaussian  noise,  at  each  of  4  known  variances  were 
generated.  These  sequences  were  then  passed  through  the  optimum  detector  for 
this  case  (signal  known  exactly  in  noise  known  exactly  (SKE  in  NKE)).  The  ROC 
curves  were  then  computed  and  drawn.  The  signal  to  noise  ratio  (SNR)  is: 


SNR  =  — 

W0 


110 


(7.9) 


where  £  is  the  signal  energy,  and  N0  is  the  noise  power  spectral  density. 

Figure  7.7  contains  the  ROC  curves  for  the  4  different  SNR’s.  These  match  well  the 
ROC  curves  originally  reported  in48  for  this  case: 
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Figure  7.7  SKE  in  Gaussian  NKE  ROC  curves 

l  wo  different  approaches  are  considered  for  the  Bayes  optimum  detector, 
namely  block  processing  and  time  sequential.  These  were  described  in  detail  in 
Chapter  6,  and  the  next  two  chapters  are  concerned  with  the  detector  comparisons 
subject  to  these  approaches.  We  now  posess  the  tools,  and  are  ready  to  examine 
the  detector  performance  evaluations  and  comparisons. 
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8.  Block  processing  optimum  detector  vs.  ad-hoc 
detector  comparisons. 

In  this  chapter,  the  performance  of  a  block  processing  Bayes  Optimum 
Detector  (BOD)  is  compared  to  that  of  an  ad-hoc  adaptive  detector.  The  reason  the 
(sub-optimum)  block  processing  approach  is  considered  here  at  all  is  that  for 
i^al-time  systems,  this  sometimes  is  the  only  viable  approach.  The  test  scenario 
under  which  both  detection  approaches  were  tested  was  described  in  Chapter  7. 
The  array  used  in  the  comparison  is  a  four  element  array,  where  the  main  beam 
was  composed  of  summing  all  four  element  outputs,  and  the  reference  beam(s) 
was  generated  by  subtracting  two  adjacent  element  outputs.  When  the  fixed 
interference  is  present,  the  problem’s  order  is  increased,  and  there  is  a  need  for 
another  degree  of  freedom.  Then,  another  reference  beam  is  generated,  and  both 
detector  schemes  tested  have  another  input  available  for  processing. 

Figure  8.1  depicts  the  typical  range-Doppler  map  of  a  synthesized  ping 
used  in  the  basic  test  scenario  where  the  interference  direction  is  assumed  known 
precisely  (i.e.  the  sonar  array  is  at  a  known  depth).  Here,  since  the  sonar  array  is 
in  relatively  shallow  water,  surface  reverberation  is  hard  to  distinguish  from 
volume  reverberation  since  they  both  have  an  almost  zero  Doppler  shift.  Figure 
8.2  shows  the  adaptive  filter’s  output  range-Doppler  map. 
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Figure  8.3  plots  main  beam  power,  reference  beam  power,  adaptive  filter's  output 
power  and  reverberation  cancellation  power  (ratio  between  adaptive  filter's  output 
power  and  main  beam  power).  Here,  one  can  see  how  much  of  the  surface 
interference  is  canceled  by  the  ad-hoc  detection  structure's  adaptive  filter  around 
0.2  sec,  where  both  detectors'  outputs  are  sampled. 


Main  been  MqrtJue  output 


CL  1  *.4  W.M  w.w  w-v 
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Figure  8.3  Main,  reference  and  filter's  output  power. 


8.1.  Known  depth 

Each  comparison  between  the  two  detection  approaches  is  based  on  a 
Monte-Carlo  simulation  of  500  runs.  In  the  first  comparison  made,  the  sonar  array 
was  placed  at  a  known  depth,  and  the  Bayes  optimum  detector  was  given  that 
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depth.  Signal  to  surface  interference  ratio  was  set  to  2.8.  Figure  8.4  shows  the  ROC 
curves  of  both  detectors. 


5KE  in  NKD  _  Boyes  Opiimum  Detector  vs. 
Adaptive  Noise  Cancel er.  Alpha  =0.02 


ROC 
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Figure  8.4  The  two  detectors’  ROC  curves  (ANC  is  *****) 

It  is  obvious  that  the  Bayes  Optimum  detector  (BOD)  performs  better  than  the  ad- 
hoc  Adaptive  Noise  Canceler  (ANC).  in  fact,  compare  these  ROC  curves  to  the 
curves  given  in  Figure  7.7  (theoretical  SKE  in  Gaussian  NKE  (SKENKE)),  which  can 
serve  as  a  reference.  One  sees  that  the  BOD  performs  close  to  a  SKENKE  detector 
when  SNR-0.5,  while  the  ANC  performs  as  a  SKENKE  detector  when  SNR-0.125. 
This  translates  into  a  "6  dB  performance  difference. 
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8.2.  Uncertain  depth 

In  the  next  comparison,  array  depth  is  not  a  known  quantity.  What  is 
known  is  its  probability  density  function,  which  is  assumed  uniform  between 
16.5m  and  43.5m.  In  this  way,  the  average  depth  is  kept  the  same  as  the  previous 
comparison’s  depth  (30m).  The  comparison  is  based  again  on  500  Monte-Carlo 
simulation  runs,  where  the  array  was  physically  placed  at  various  depths, 
commensurate  with  a  uniform  distribution.  Figure  8.5  shows  the  ROC  curves  for 
this  comparison. 


SKE  in  NUO  -  27  Meiers  Uncertainty 
Bayes  Optimum  Detector  vs. 
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Figure  8.5  Uncertain  depth  (27m)  ROC  curves  (ANC  is  *+***) 

It  is  clear  that  the  ANC’s  performance  hasn’t  changed,  since  it  has  the  ability  to 
adapt  to  the  varying  depth.  The  BOD's  performance  has  degraded,  but  with  this 
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level  of  uncertainty,  it  is  still  performing  better  than  the  ANC.  Naturaly,  we  will 
now  check  what  happens  when  the  depth  uncertainty  is  increased.  Depth  was 
given  twice  the  uncertainty  (i.e  54m),  still  keeping  the  average  depth  at  30m  as 
before.  Figure  8.6  plots  the  ROC  curves  for  this  case. 
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Figure  8.6  Uncertain  depth  (54m)  ROC  curves  (ANC  is  *****) 

Here,  the  BOD  has  basically  fallen  apart,  since  the  interference  direction 
uncertainty  is  large.  The  ANC,  whose  performance  is  essentially  unchanged,  is  now 
performing  better  than  the  BOD,  using  its  capability  to  adapt  to  the  varying 
interference  direction. 


Another  interesting  point  is  to  see  what  happens  when  the  signal  to  noise 
ratio  is  increased.  This  was  done  for  the  27m  of  depth  uncertainty,  and  an  SNR  of 
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11.2.  Figure  8.7  shows  the  comparison  result. 
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Figure  8.7  Higher  SNR  ROC  curves  (ANC  is  *****) 

As  can  be  seen,  the  performance  of  both  detectors  has  improved  as  expected. 

8.3.  Uncertain  depth  plus  a  fixed  interference 

At  this  point  in  time,  it  was  clear  that  as  soon  as  the  depth  (interference 
direction)  uncertainty  exceeds  some  value,  the  ANC  starts  performing  better  than 
the  BOD.  The  next  step  was  to  check  what  effect  an  extra  interference  source  has 
on  both  detectors'  performance.  An  interference  source  whose  direction  is 
uncertain,  but  whose  direction  probability  density  function  was  known,  has  been 
added  into  the  scenario.  The  interference  was  coming  from  directions  between 


1 19 


-10"  and  -30“  .  Interference  power  level  was  made  equal  to  the  boundary 
interference  power  at  the  detection  range.  Figure  8.8  shows  the  results,  we  see 
that  the  result  of  the  added  interference  is  that  the  BOD  performs  poorly. 
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Figure  8.8  Fixed  interference  introduced  (ANC  is  +****) 
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8.4.  Conclusion 

When  comparing  Figure  8.8  to  Figure  8.5  (see  Figure  8.9),  we  see  that  the 
ANC's  performance  remained  essentially  the  same,  even  after  the  addition  of  the 
fixed  interference.  This  demonstrates  its  ability  to  adapt  to  another  interference 
source,  when  the  adaptive  filter  order  is  given  another  degree  of  freedom.  As 
expected,  the  BOD  performs  very  poor,  and  the  effect  another  interference  source 
has  on  it,  is  akin  to  the  effect  that  an  increased  depth  uncertainty  had  (see  Figure 
8.6).  Figure  8.9  combines  Figures  8.5  and  8.8,  to  compare  both  detectors’ 
performance  under  the  same  depth  uncertainty  (27m),  with  and  without  the  fixed 
interference. 
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Figure  8.9  Performance  with  and  without  fixed  interference  (ANC  is  ♦****) 

We  have  compared  the  performance  of  the  two  detection  approaches, 
when  the  Bayes  Detector  was  restricted  to  block  processing.  Since  this  restriction 
makes  the  detector  sub-optimal,  its  performance  has  not  been  uniformly  better 
than  its  ad-hoc  counterpart.  The  next  chapter  will  remedy  that  by  using  a  time 
sequential  approach  for  the  BOD. 


9.  Time  sequential  optimum  detector  vs.  ad-hoc 
detector  comparisons. 


In  this  chapter  we  will  remove  the  block  processing  restriction  from  the 
Bayes  Optimum  Detector  (BOD)  and  let  it  make  use  of  all  the  available  data.  This 
way,  the  BOD  becomes  a  truly  optimum  solution,  and  the  comparison  between  the 
two  detectors  is  more  fair.  This  will  also  enable  us  to  demonstrate  the  adaptive 
nature  of  the  BOD.  The  difference  between  the  block  processing  and  the  time 
sequential  approaches  has  been  detailed  in  Chapter  6. 

9.1-  Uncertain  depth 

The  first  comparison  is  identical  to  the  second  comparison  made  in 
Chapter  8  (SKE  in  NUD).  Here,  there  is  a  27m  depth  uncertainty,  and  the  BOD  is 
time  sequential.  Figure  9.1  shows  the  ROC  curves  for  this  case. 
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Figure  9.1  Time  sequential  BOD  vs.  ANC  (ANC  is  *+***) 

Obviously,  the  performance  difference  between  the  BOD  and  the  ANC  has 
increased  (compare  to  Figure  8.5).  Now  that  the  BOD  is  making  use  of  past  data  to 
update  the  a-priori  pdf,  it  performs  better.  As  before,  the  next  step  is  to  increase 
the  depth  (interference  direction)  uncertainty,  and  observe  the  performance.  The 
uncertanty  is  increased  to  54m  as  before,  and  Figure  9.2  details  the  result. 
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Figure  9.2  54m  uncertainty  ROC  curves  (ANC  is  **♦**) 

As  opposed  to  what  happened  in  the  same  situation  when  the  block  processing 
approach  was  used,  now  even  with  the  increased  depth  uncertainty,  the  BOD  still 
performs  better  than  the  ANC,  though  the  performance  differential  is  reduced. 

As  described  in  Chapter  6,  the  time  sequential  detector  presents  an 
adaptive  quality  in  that  it  updates  the  a-priori  probability  density  functions  (pdfs) 
of  the  uncertain  parameters  from  iteration  to  iteration,  thereby  carrying  important 
information  from  the  past  to  the  present.  It  is  interesting  to  watch  how  these  pdfs 
adapt  as  time  progresses  (see  equations  (6.8)  and  (6.9)).  Figure  9.3  depicts  the 
behaviour  of  the  array  element  to  element  delay  (AEED)  pdf  for  the  iterations 
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preceding  the  one  in  which  detection  is  made.  AEED  is  a  linear  function  of 
interference  direction,  and  was  the  unknown  parameter  used. 
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Updating  af  the  a-priori 
probability  density  Function 


Bayes  Optimum  Time  Sequential  Detector 


Figure  9.3  Updating  pdf’s. 
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Note  how  system's  knowledge  of  the  AEED  improves  from  very  diffuse  (Iteration 
#1),  to  more  and  more  precise,  as  time  progresses.  Since  the  surface  return  comes 
from  shallower  and  shallower  angles  as  the  ping  propagates  through  the  water,  the 
AEED  becomes  smaller  and  smaller,  and  so  does  its  estimate. 

9.2.  Uncertain  depth  plus  a  fixed  interference 

As  before,  we  now  check  what  effect  an  extra  interference  source  has  on 
both  detectors'  performance.  An  interference  source  whose  direction  is  uncertain, 
but  whose  direction  probability  density  function  is  known  has  been  added  into 
the  scenario.  The  interference  is  coming  from  directions  between  -10°  and  -30°. 
Interference  power  level  was  made  equal  to  the  boundary  interference  power  at 
the  detection  range.  Figure  9.4  shows  the  results. 
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SKE  fn  NUD  -  27  Meiers  Uncertainty 
plus  a  Fixed  Interference 
Bayes  Optimum  Time  Sequential  Detector 
vs.  Adaptive  Noise  Cancel er 


ROC 


Figure  9.4  Fixed  interference  introduced  (ANC  is  *****) 

Here,  the  BOD’s  performance  has  degraded  due  to  the  presence  of  the  fixed 
interference.  Even  though  the  ANC  seems  to  be  less  sensitive  to  the  presence  of 
the  fixed  interference,  the  BOD  still  performs  better  than  the  ANC.  This  is 
contrasted  with  Figure  8.8,  where  the  ANC  performs  better  than  the  block 
processing  BOD. 


9.3.  No  surface  interference 

It  is  interesting  to  check  what  influence  removing  the  surface  interference 
from  the  scenario  will  have  on  both  detectors’  performance.  Time-varying 
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boundary  interference  power  and  fixed  interference  power  were  purposefully  kept 
at  an  equal  level.  It  is  expected,  therefore,  that  performance  will  be  very  close  to 
the  case  where  only  the  surface  interference  was  present.  Figure  9.5  shows  the 
results: 


5KE  in  NUD  -  27  Meiers  Uncertainty 
plus  a  Fixed  Interference 
no  surface  reverberation 
Bayes  Optimum  Time  Sequential  Detector 
vs.  Adaptive  Noise  Cancel er 
ROC 


pr 

Figure  9.5  Surface  interference  removed  (ANC  is  *****) 

It  is  clear  that  indeed  the  performance  is  very  similar  to  the  one  reported  in  Figure 


Lastly,  Figure  9.6  shows  the  detector  performance  when  depth  uncertainty 
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is  increased  from  27m  to  54m,  and  when  the  fixed  interference  is  present. 


5KE  in  NUD  -  54  Meiers  Uncertainty 
plus  a  Fixed  Interference 
Bayes  Optimum  Time  Sequential  Detector 
vs.  Adaptive  Noise  Cancel er 


ROC 


pr 

Figure  9.6  Higher  depth  uncertainty  (ANC  is  *****) 

The  performance  of  both  processors  has  degraded  uniformly,  while  the  BOD  still 
performs  best,  as  expected. 


9.4.  Conclusion 

When  a  truly  optimum  approach  is  taken,  and  the  Bayes  optimum  detector 
is  allowed  to  make  use  of  all  the  available  data,  it  performs  better  than  the  ad-hoc 
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adaptive  noise  canceler.  With  the  time-sequential  approach,  the  BOD  also  posesses 
an  adaptive  quality,  in  updating  the  a-pricri  pdf’s  of  the  unknown  parameter  from 
iteration  to  iteration. 


Conclusion 

The  active  sonar  problem  has  been  carefully  examined  in  this  dissertation. 
Specifically,  situations  under  which  the  sonar  suffers  from  boundary  interference 
were  analyzed.  Underwater  acoustic  reverberation,  which  is  typically  the  limiting 
noise  process  for  this  problem,  has  been  studied.  Then,  two  main  schemes 
designed  to  detect  a  desired  signal  buried  in  ocean  reverberation  were  suggested 
and  examined.  These  two  schemes  represent  two  different  approaches  to  the 
problem.  The  first  is  an  ad-hoc  engineering  approach,  which  uses  an  adaptive 
noise  canceler  in  the  front  end,  to  get  rid  of  the  undesired  boundary  interference, 
and  then  uses  a  matched  filter  assuming  now  that  the  problem  is  one  of  known 
Gaussian  statistics.  The  second  approach  is  a  more  global  one.  It  makes  use  of  the 
known  statistical  properties  of  the  problem  to  arrive  at  an  optimal  solution  in  the 
Bayesian  sense. 

Recognizing  that  close  analyse  expressions  for  system  performance  are  often 
very  hard  to  derive,  an  important  tool  has  been  developed.  This  software  package 
is  capable  of  synthesizing  multi-element  (or  beam)  time  series  outputs,  taking  into 
account  transmit  and  receive  beampatterns,  the  ocean  environment  including 
surface,  volume  and  bottom  scattering  strengths,  and  the  sonar  platform 
dynamics.  It  sets  up  the  normal  equations  for  the  problem,  and  uses  the  solution 
to  create  a  time  varying  HR  filter  which  represents  the  time  varying  characteristics 
of  a  sonar  ping  propagating  in  the  ocean.  Once  the  filter  has  been  derived,  an 
ensemble  of  sonar  pings  may  be  collected  by  passing  uncorrelated  white  noise 
vectors  through  the  filter.  This  tool  is  then  heavily  used  in  assessing  performance 
of  the  various  detection  algorithms. 

Next,  a  multi-channel,  joint  process  pole-zero  adaptive  filter  is  developed  and 
its  performance  is  compared  to  its  all-zero  counterpart.  It  is  shown  that  some 
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cases  exist  where  the  ability  to  adapt  poles  is  important. 

Last,  but  not  least,  the  performance  of  the  ad-hoc  detectors  is  compared  to 
that  of  Bayes  optimum  detectors.  It  is  shown  that  optimum  detectors,  when 
allowed  to  make  use  of  all  the  available  data,  perform  better  than  the  ad-hoc 
adaptive  detectors.  Moreover,  it  is  shown  that  the  optimum  detector  possesses  an 
adaptive  learning  quality,  through  the  process  of  updating  the  a-priori  probability 
density  functions  of  the  unknown  parameters. 

The  work  done  in  this  dissertation  can  be  extended  in  a  few  directions.  The 
multi-channel  joint  process  ARMA  adaptive  filter  developed  here  can  be  used  in 
various  noise  canceling  and  spectral  estimation  problems,  and  its  ability  to  adapt 
poles  can  be  exploited.  The  reverberation  simulation  tool  (REVSIM)  incorporates 
only  isovelocity  sound  speed  profiles,  and  the  next  step  would  be  to  allow  a 
variable  sound  speed.  Other  reverberation  synthesis  techniques  exist,  and  have 
been  implemented  by  other  researchers.  It  would  be  interesting,  and  beneficial  for 
the  active  sonar  community,  to  compare  and  evaluate  the  different  methods 
against  each  other.  The  application  of  Bayes  optimum  detectors  when  there  are 
uncertain  parameters  has  been  outlined  and  demonstrated  for  the  active  sonar 
problem.  Two  uncertain  parameters  concentrated  on  were  platform  depth,  and 
direction  of  arrival  of  a  fixed  interference.  It  would  be  very  interesting  to  extend 
the  methods  and  ideas  outlined,  to  other  uncertain  parameters  such  as  sound 
velocity.  Finally,  the  investigation  of  applying  Bayes  optimum  schemes  for  feature 
estimation  and  extraction  (e.g  direction  of  arrival)  using  their  adaptive  learning 
capability,  may  prove  beneficial. 

i 

I 
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Appendix  A 


The  linear  minimum  mean  square  estimation 
problem: 

The  problem  is  presented  in  the  following  figure: 


Figure  A.l  -  The  MMSE  estimation  problem 

a(n)  is  a  causal  FIR  filter  of  length  p+\.  The  goal  is  to  filter  the  stochastic 
time  series  x(n)  to  yield  an  optimal  estimate  (in  some  sense)  of  the  "desired"  time 
series  d(n).  The  error  is  given  by: 

p 

ed(n)  =  d(n)+Y^akx(n-k)  (A.l) 

o 

or  in  vector  notation: 


ed(n)  =  d(n)+aT x  (A.2) 

where  r  represents  the  vector  transpose  operator.  The  solution  to  the  problem 
depends  on  the  assumptions  made  about  the  input  time  series,  and  the  optimality 
criterion  chosen.  The  Wiener  filtering  approach  assumes  the  input  time  series  is 
|  stochastic,  and  seeks  to  minimize  the  mean  square  error: 

i 

) 

min£ 


ed(n) 


(A.3) 
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where  £[  .  ]  is  the  expected  value  operator. 


eJ(n)  |  2  =  |  d(n)-t-aTx  i  2 

=  d(n)d  {n)+d(n)a‘ x.+an x  {n)d(n)+an  x  (n)a'x(n) 


(A. *4) 


where  H  represents  the  Hermitian  operator.  Now,  taking  the  expected  value  yields: 


£  I 


ed(n) 


t*(a)  =  0|^)+^Hi2+i2HJ+i2H[4>]ri2 


(A.5) 


where  od0  is  the  desired  signal  power,  g  is  the  crosscorrelation  vector  between 
the  input  x(n)  and  the  desired  signal  d(n),  and  [<t>]  is  the  autocorrelation  matrix  of 
the  input  process  x(n).  Minimizing  E^ia)  with  respect  to  a  yields: 


or 


0  =  a+[<b]Ta 


(A.6) 


[<t>]T<z  =  -a 

which  are  called  the  normal  equations.  Solving  for  a  we  get: 


(A.7) 


a 


(A.8) 


In  the  case  that  the  input  time  series  x(«)  is  stationary,  [<t>]  is  of  a  Toeplitz  form: 


=  tj-k  (A-9) 

Here,  all  matrix  elements  along  any  diagonal  are  equal.  This  becomes  useful  when 

solving  (A.7),  avoiding  the  need  to  directly  invert  [<f>], 
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Appendix  B 


Transition  equations  for  the  time-varying  element  to 
element  delay. 

In  time-sequential  processing,  when  transitioning  from  one  iteration  to  the 
next,  an  update  is  made  on  the  a-priori  pdf’s.  Since  <5,,  the  element  to  element  time 
delay  is  time-varying,  a  transition  equation  is  required,  which  determines  the  next 
iteration's  time  delay  as  a  function  of  this  iteration's  delay.  Consider  the 
problem's  geometry  as  is  shown  in  figure  B.l: 


Ocean  surface 


Figure  B.l  -  element  to  element  time  delay  computation 


At  time  T, ,  the  element  to  element  delay  is: 


m 


Id 

cZ(Tj) 


Id 


_K 

T, 


(B.l) 
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where 


At  time  T(+1,  the  element  to  element  delay  is: 


5,+i(7i+1) 


K 

Tm 


(B.2) 


where 


T;+1  -  Tj+AT 

Let 


Si+l  = 


(B.3) 


then 


A  =  6l+l-S,  = 


K  KTj—KTt—K&T 


7,+A7  Tj  7,(7,+ A  7) 

Substituting  A  back  into  (B.3)  yields: 

5/+l  = 


KA7  AT 

- --V»+i -  (B.4) 

T,(Tj+ A  71  AC 


AT 

AC 


(B.5) 


«/+1  ~ 


A  7 

1+5, - 

K 


Now,  since 


d 

T-,  =  -+(/~l)A7 
c 


then  (B.5)  becomes: 
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Now, 


*,+.= 


«/ 


5/ 


1+- 


K  AT 


T,  K 


AT 

1  + - 

T, 


1+- 


A  T 


'+(/— 1)AT 


5/  = 


Id 


Id 


cT, 


+(/—  1)AT 
1  c 


d  = 


cZ6,(i-l)AT 

l—cs, 


;  i>l 


(B.6) 


(B.7) 


Substituting  (B.7)  back  into  (B.6)  yields: 


6, 


6i+ i  - 


*/ 


1+- 


1 


C5f(/ —  1) 

/-C«, 


-+(/-!) 


l—CSj 

id- 1) 


i+- 


;  />i 


(B.8) 


which  is  the  desired  result  used  in  (6.10). 


Part  of  the  update  equation  for  the  a-priori  probability  density  function  of 


6,  is  the  derivative 


d5l+]  d8i 

- I  5  .  It  may  be  more  convenient  to  compute  |  - 

ds,  1  d&M 


and  then  use: 


di, 


i+i 


ds,  d5j 


d5; 


+1 


Using  (B.8)  we  have: 


139 


Si 


~  ^/+I  +  " 


i+l 


ils, 


i+l 


C5((i-l) 

l-CSf 


'+(1—1) 


/(/-1)+C5,+1 


d6j  //(/(/- l)+c5f^)-ri/5,+1  il2(i-l) 

1  (/(/— l)+C<5;+i)  (/(/—  l)+c5,+1) 

Now,  inverting  the  derivative  we  get: 


1  (1(1—1  )+c^+[) 

d*i  //2(/-l) 

and  finally: 


c/5 


/(/-!)+- 


C5,. 


1+' 


i+i 


l-cg/ 

/(l-l>  ' 


c/5, 


»/2(/-l) 


(B.9) 


(B.IO) 


(B.l  1) 


(B.12) 
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